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ABSTRACT 

' We present a method that self-consistently tracks the growth of supermassive black 

I holes (BHs) and the feedback from active galactic nuclei (AGN) in cosmological, hy- 

• , drodynamical simulations. Our mo del is a substantially modified version of the one 

' introduced by ISpringel et al.l (|2005[ ) implemented in a significantly expanded version 

of the GADGET III code, which contains new prescriptions for star formation, super- 



Q ■ nova feedback, radiative cooling and chemodynamics. We simulate the growth of BHs 

^ ' from an initial seed state via Eddington-limited accretion of the surrounding gas, and 

C/3 , via mergers with other BHs. Because cosmological simulations at present lack both 

^ , ' the resolution and the physics to model the multiphase interstellar medium, they tend 

to strongly underestimate the Bondi-Hoyle accretion rate. To allow low-mass BHs to 
. grow, it is therefore necessary to increase the predicted Bondi-Hoyle rates in star- 

^ ' forming gas by large factors, either by explicitly multiplying the accretion rate by a 

I numerical correction factor, or using an unresolved, subgrid model for the gas close to 

i the BH. We explore the physical regimes where the use of such multiplicative factors 

' is reasonable, and through this introduce a new prescription for gas accretion by BHs. 

Feedback from AGN is modeled by coupling a fraction of the rest-mass energy of the 

■ accreted gas thermally into the surrounding medium. We describe the implementation 
' as well as the limitations of the model in detail and motivate all the changes relative 

0^ . to previous work. We demonstrate how general physical considerations can be used to 

■ choose many of the parameters of the model and demonstrate that the fiducial model 

reproduces observational constraints. 

We employ a large suite of cosmological simulations, in which the parameters of 
' the BH model are varied away from their fiducial values, to investigate the robustness 

M I of the predictions for the cosmic star formation history and the redshift zero cosmic 

BH density, BH scaling relations, and galaxy specific star formation rates. We find 
that the freedom introduced by the need to increase the predicted accretion rates 
by hand, the standard procedure in the literature, is the most significant source of 
uncertainty. Our simulations demonstrate that supermassive BHs are able to regulate 
their growth by releasing a fixed amount of energy for a given halo mass, independent 
of the assumed efficiency of AGN feedback, which sets the normalization of the BH 
scaling relations. Regardless of whether BH seeds are initially placed above or below 
the BH scaling relations, they grow onto the same scaling relations. AGN feedback 
efficiently suppresses star formation in high-mass galaxies. 
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Over the past decades a growing body of observa- 
tional and theoretical evidence has suggested that su- 
permassive black holes (SMBHs; niBH > lO^M©) ex- 
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ist in the centers of all ga l axies with spheroid s (e 
iKormendv fc Richston3 Il995l : iFerrarese fc MerrittI [20oi 
and that the properties of these SMBHs are tightly cor- 
related with the properties of the spheroid in which 
they reside. For example, the mass of the SMBH is 
found to be tig htly correlated w i th th e bulge stellar mass 
or luminosity (Maeorrian et al. Il998l : [ McLure fc Dunlod 
l2002l : iHaring fc Rix 2004; 'Laor 200 ll), stellar velocity dis- 
persion jGebhardt et al, ,2000, ; Merritt fc Ferraresd I2OOII ; 
iTremaine et alj 12002! ) , ajid galaxy concentration, a s mea- 
sured by the Sersic index (|Graham fc Driverll2007l '). Some 
recent work has demonstrated that these correlations may 
be understood in terms of a black hole (BH) 'funda- 
mental plane', relating BH mass, galaxy effective radius, 
stellar velocity dispersion and stellar mass. Here, the 
mass of the SMBH essentially tracks the bindin g energy 
of th e stellar bulge |Ma rconi fc Hunt 2003; Fcol i fc Melel 
I2OO5I ; lAUer fc Richstonel [20071 : iHopkins et all I20"0^ ). al- 
though other authors argue that the appearance of a funda- 
mental plane is actually d ue to biasing c aused by the pres- 
ence of galaxies with bars l|Grahamll200i ). 

The exact mechanisms leading to the tight observed 
coupling between galaxy spheroidal components and cen- 
tral active galactic nuclei (AGN) are not yet fully under- 
stood, although it has long been r ecognized that the for- 
mation m echanisms of SMBHs (e.g. ISilk fc Ree s' 1998) and 
stars (e.g. lDekel fc Sil5il986 ) are most likely self-regulating. 
These results suggest that the same processes that shape 
galaxy spheroids also act on the central BHs. Correlations 
between AGN activity and other processes provide other 
clues about the mechanisms that lead to the buildup of the 
SMBH population. There is evidence that there exists a link 
between galactic star formation and accretion onto a cen- 
tral AGN: in a global sense , t he e volution of the cosmic star 
formation rate (e.g. ,Ma dau et a l.; 1 9961 ) and the luminosity 
densi ty of quasars are tightly correlated l|Bovle fc Terlevicbl 
Il998l ). Additionally, on the scale of individual objects it has 
been found that the most powerful narrow line AGN are 
preferentially found in galaxies that appear to have under- 
gone a recent starburst phase (Kauffmann ct al. 200j). 

The massive BHs present in the centres of galaxies are 
likely to have started their lives as 'seed' BHs. The typi- 
cal masses of seed BHs remains somewhat uncertain, and 
depends upon the mechanism by which they form. Plausi- 
ble mechanisms include the collapse of population HI stars, 
giving rise to BHs with masses in the range 10'^ Mp < 
rwRH < IO^Mm (e.g.lMadau fc Reesll200ll : llslam et al.ll2003l ; 
[Schneider et al.ll200a r and direct collapse of matter in high 
redshift, low angular momentum haloes, whic h may give 
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rise to seed BHs with masses ~ 10^ Mp (e.g. Loebl 19941: 



Bromm fc Loebll2003l: iBegelman et al.ll2006l : iDiikstra et al] 
20081 : IVolonteri fc Nataraianll2009l ). These seed mass BHs 
can then grow either by mergers with other BHs (e.g. 
llslam et alTl2003l ). or through accretion of gas and/or stars. 

Accretion of matter onto central BHs, accompanied by 
the release of a fraction of the rest mass energy of the fuel, 
has long been recogniz ed as one of t he most likely mech- 
anisms to power AGN l|Salpeteij|l964l ). and a coupling be- 
tween the accretion history of an AGN and the gas dynamics 
of the bulge provides a plausible mechanism by which AGN 
and bulge properties could become strongly correlated cou- 
pled (e.g. lSilk fc Ree3ll998h . For example, it has been sug- 



gested that the central BHs grow until they release suffi- 
cient energy to unbind t he gas that feeds them from the 
host galaxy (Fabian| |l999l ). Bursts of AGN activity then ex- 
pel gas from galaxies and remain quiescen t until stellar mass 
loss r eplenishes the galaxy's gas reservoir (jCiotti fc Ostrikerl 
I2OOII) 

A theoretical link between galaxy mergers and both 
galaxy-scale starburst events and active AGN phases 
has been well established and modelled. Galaxy merg- 
ers have long been recognized as a mechanism by which 
gas ca n potentially be chann eled to the centre of a 
galaxy (|Toomre fc Toomrelll97j ). and N-body simulations 
of galaxy mergers confirmed and extended this picture 
by showing that the asymmetrical gravitational potential 
present during mergers is capable of funneling gas ef ficiently 
to the center of a galaxy l|Mihos fc Hernauistlll994l ). where 
it may be accreted by a SMBH. Hydrodynamical s i mula- 
tions of galaxy mergers (iBarnes fc Hernguistl Il99ll . Il996l : 



iMihos fc Hernauist|[l996l : iKapferer et al.ll20q5l), an d numer- 
ical models of AGN growth (|Kapferer et al T I2OO5I : ISprin"geil 
l2005h predict that these merger events are indeed responsi- 
ble f or the rapid grow th of AGN. Recent numerical studies 
(e.g. iMicic et al.ll2007l ) have indicated that both BH merg- 
ers and gas accretion are important processes in forming the 
population of BHs that we observe in the local universe. 

Thus, it seems that we can paint a coherent picture 
in which emission by AGN, galaxy mergers and the growth 
of supermassive BHs are closely intertwined. As such the 
study of any one of these processes requires an under- 
standing of all of them. For this reason detailed studies 
of the co-evolution of the AGN and galaxy populations 
in a cosmological context often resort to numerical tech- 
niques. Early theoretical studies of the galaxy-AGN con- 
nection relied upon dark matter halo merger rates, with- 
out any separate treatme nt of galaxy forma tion p rocesses 
lEfstathiou fc Reed 1 19881 : iHaehnelt fc Reed Il993l ). Later 
work expanded upon this groundwork by incorporating 



mation (e.g. 


Kauffmann fc Haehnelt' '2OOO'; 'Cattanco' 2001 


Benson et al. 


[2003: Granato et al. 2004; Croton ct al. 2006 


Bower et al.l 


20061: iLagos et al.l I2OO8I: iBaughl I2OO6I). Semi- 



analytic models indicate that feedback from AGN is nec- 
essary in order to build up a red-sequence of galaxies. Al- 
though a combination of photo-heating by reionization and 
supernova feedback can suppress star formation in low mass 
haloes - bringing the galaxy luminosity function in line with 
observations at the low mass end - models without AGN 
feedback face the problem that the reheated gas in mas- 
sive haloes would eventually coo l, giving rise to an excessive 
number of bright galaxies (e.g. iBower et al.l 12009 ) as com- 
pared to the local universe. 

A more computationally challenging approach is to sim- 
ulate galaxies hydrodynamically, with additional sub-grid 
modelling of the growth and energy feedback from AGN. 
Numerical hydrodyn amic simulations of galaxy mergers 
conta i ning AGN (e.g. ISpririgel et al.ll2005l: iDi Matteo et all 
I2OO5I ; iHopkins et all l2006l : [Robertson et al.l |20 06 l) haw 
shown that the presence of a central AGN can significantly 
alter the structure of merger remnants, particularly by ex- 
pelling a hot halo of diffuse, low-angular momentum gas 
from the center of the remnant. More recent numerical stud- 
ies have revealed that dissipation and dry mergers are likely 
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to play a fundam ental role in shapin g the co-evolution of 
BHs and galaxies (|Hopk ins et al ] |200^ Hydrodynamic sim- 
ulations of fuU cosmological vol umes (Di Mattco ot al. '2008"; 
ISiiacki et al] l2007l : ICroft et al.| [2008: Okamoto et al. 2008) 
have probed the effect of AGN on a cosmologically represen- 
tative set of galaxies and showed that the inclusion of AGN 
physics into galaxy formation simulations allows us to match 
many of the observed properties of galaxies in the local uni- 
verse. The modelling of an AGN population in this manner 
is both computationally very expensive and subject to very 
many, as yet poorly understood, numerical effects. As such 
studies of this type must take care to test the robustness of 
the models to all physical and numerical parameters. 

The focus of the current work is to present and test a 
new model for the co-evolution of BHs and galaxies. We note 
that nearly all BH models published thus far in the literature 
e mploy the star formation and supernova feedback models 
of ISpringel fc Hernguisj l|2003l) (hereafter S H03), and the 
model for BH growth and AGN feedback of ISpringel et al] 
120051 ^1 (hereafter S05). Throughout this paper we high- 
light similarities and differences between our approach and 
that used previously in the literature. The primary differ- 
ence between our models and previous work is that we 
employ a different parametrization of the process of gas 
accretion onto BHs as well a different implementation of 
AGN feedback. We show that changes to the BH accretion 
model can lead to profound differences in galaxy proper- 
ties, global star formation rates and BH demographics. We 
examine both the global properties of the simulation, such 
as the integrated star formation rate and cosmic BH den- 
sity, and consider the properties of individual galaxies, in- 
cluding specific star formation rates, and the BH funda- 
mental plane. We quantify how uncertainties in our nu- 
merical model and all of our parameter choices affect the 
reliability of our results. We find that changes in the nu- 
merical model that generates seed mass BHs, and in the 
model that distributes feedback energy into the ISM do not 
strongly affect our results. However, the accretion model 
is found to be of crucial importance in understanding our 
results. Throughout we assume a flat ACDM cosmology 
with the cosmological parameters: {f2,„, fib, Q.h, erg, n^, ft} = 
{0.238, 0.0418, 0.762, 0.7 4,0.951,0.73} as d etermined firom 
the WMAP 3-year data (ISpergel et al.ll2007l'l and co nsistent 
with the WMAP 5-vear data l|Komatsu et al.ll2008h . Where 
necessary, observational results have been scaled to our cho- 
sen cosmology, and the stellar initial mass function (IMF) 
assumed in observational analyses has been scaled to the 
Chabrier IMF used in our simulations. 

The paper is structured as follows: In Sec. [5] we intro- 
duce our simulation set, and describe briefly the sub-grid 
physics modules that are not directly related to BHs. In 
Sec. 13] we describe in detail our model for BH formation, 
growth and feedback and we motivate our choices for nu- 
merical parameters in Sec.|4l In Sec.[5]we present simulation 
results, including a comparison with redshift zero observa- 
tional data and an investigation into the severity of uncer- 
tainties introduced by different parameter choices. Finally, 



Although see lOkamoto et al. for a different approach. 

^ Our value of trg is 1.6a lower than allowed by the WMAP 5-year 
data. 



in Sec. |6] we discuss and summarize our findings. In a com- 
panion work we investigate in detail the interplay between 
feedback from AGN and other feedback processes, including 
winds driven by Type 2 supernovae and mass loss from the 
stellar population. 



2 NUMERICAL SIMULATIONS 

In this section we introduce the numerical techniques used in 
our simulations and provide a brief overview of the sub-grid 
physics modules that are not directly related to BH growth 
or AGN feedback. 

We have carried out a suite of cosmological simulations 
using Smoothed-Particle Hydrodynamics (SPH ) ||LucvI1977| : 
iGingold fc Monaghanlll977j : lMonaghanlll992l '). employing a 
significantly exte nded version of the parallel PMTree-SPH 
code GADGET III ('Springcl"2005':'Sprin gel. Yoshida fc White! 
l200ll ). a Lagrangian code used to calculate gravitational and 
hydrodynamic forces on a particle- by-particle basis. The 
initial particle positions and velocities are set at z = 127 
using the Zel'dovich approximation to linearly evolve posi- 
tions from an initially glass-like state. The production sim- 
ulations used in this study are run in boxes of size 50 co- 
moving Mpc//i, and contain 256^ particles of both gas and 
dark matter. Comoving gravitational softenings are set to 
1/25 of the mean comoving inter-particle separation down 
to 2: = 2.91, below which we switch to a fixed proper scale 
of 2 kpc/h. The production simulations have gas particle 
masses of 8.64 X lO^Mo/ft. The boxes are evolved all the 
way to redshift zero. 

In addition to hydrodynamic forces we treat star forma- 
tion, supernova feedback, radiative cooling, chemodynamics, 
black hole accretion, and AGN energy feedback in these sim- 
ulations. 

Star formatio n is tracked in the s imulations following 
the prescription of lSchave fc Dalla Vecchia (2008). Gas with 
densities exceeding a critical density for the onset of the 
thermo-gravitational instability (hydrogen number densities 
nu ~ 10~^ — 10~^ cm~^) is expected to be multiphase and to 
form stars (|Schavell2004l ). Because we lack both the physics 
and the resolution to model the cold interstellar gas phase, 
we impose an effective equation of state (EOS) with pres- 
sure P oc p"*""" for densities nn > where rig = 0.1 cm~'^, 
normalised to P/k = 10^ cm~^K at the threshold. We use 
7off = 4/3 for which both the Jeans mass and the ratio of the 
Jeans length to the SPH kernel are independent of the den- 
sity, thus preventing sp urious fragmentation due to a lack 
of numerical resolution fechave fc Dalla V ecchiall2008l ). As 
described in :Schaye_&_Dalla Vecchia (2008), gas on the ef- 
fective EOS is allowed to form stars at a pressure-dependent 
rate that reprod uces the observed Kennicutt-Schmidt law 
(lKennicut3l 19981 ) by construction, renormalised by a factoiQ 
1/1.65 to account for the fact that it assumes a Salpeter 
IMF whereas we use a Chabrier IMF. 

Energy injection due to supernovae is included 
through kinetic feedback. We employ the prescription of 

^ This normalization factor is calculated from the asymptotic ra- 
tio of the numbers of ionizing photons pre dicted from models 
of ste llar populations with a constant SFR llBruzual &: Charloj 
I2OO3I) . 
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iDalla Vecchia fc Schavd (|2008l ). which is a variation of the 
SH03 recipe for kinetic feedback. In this prescription core- 
collapse supernovae locally inject kinetic energy and kick 
gas particles into winds. The feedback is specified by two 
parameters: Firstly, the initial mass-loading rj — rhw/m*, 
which describes the initial amount of gas put into the wind, 
TTi^ , as a function of the local SFR, m,, and secondly the 
wind velocity, We use rj — 2 and — 600 km/s, which 
corresponds to 40% of the total amount of supernova energy. 
In contrast with the models of SHQ3, the kinetic energy is 
injected locally to every star formation event and wind parti- 
cles are not temporarily decoupled from the hydrodynamics 

when they are put i nto the wind. 

As described in |w lersma et al.l l|2009bl ). we follow the 
timed release of 11 different elements from massive stars 
(Type II supernovae and stellar winds) and intermediate 
mass stars (Type la supernovae and asymptotic giant branch 
stars), assuming a Chabrier initial mass function spanning 
the range 0.1 to 100 Mp). Rad i ative c ooling was imple- 
mented following IWiersma et al.l (|2009al ff. In brief, net ra- 
diative cooling rates are computed element-by-element in 
the presence of t h e cos mic microwave background and a 
iHaardt fc Madaul (|200ll ') model for the UV /X-ray back- 
ground radiation from quasars and galaxies. The contribu- 
tions of the eleven elements hydrogen, helium, carbon, nitro- 
gen, oxygen, neon, magnesium, silicon, sulphur, calcium, and 
iron are interpolated as a function of density, temperature, 
and redshift from tables that have been precomputed using 
the publicly avai lable photo-ionizatio n package CLOUDY, 
last described bv lFerland et al.l (|l998t ). assuming the gas to 
be optically thin and in (photo-)ionization equilibrium. 



3 THE BLACK HOLE MODEL 

We now provide a detailed description of our models for BH 
formation and accretion (Sec. [33J, BH mergers (Sec. [33]) 
and energy feedback from AGN (Sec. 13.3) . Throughout this 
section we highlight and justify the aspects of our model 
that differ from previous works. We also introduce all the 
relevant parameters. In Sec. [4] we motivate our choices for 
these parameters. 



3.1 Black hole formation and accretion 

Plausible BH seed formation mechanisms lead to the cre- 
ation of BHs with masses in the range 10 -lO'^ M0, whereas 
SMBHs in the local Universe have masses of up to 10^ M0 
(see Sec [1] for a discussion). To understand the origin of the 
redshift zero BH population we therefore need to model how 
BHs can grow to the sizes observed in present-day galax- 
ies. Over the past decades a picture has emerged in which 
SMBHs are embedded in dense stellar systems in the centres 
of galaxies and increase their masses pr imarily by the accre- 
tion of gas (e.E.lBegelman fc Rees|[T978l ) . BHs may also grow 
by mergers with other BHs, or by t he disruption and cap- 
ture of stars (e.g. iLvnden-Belll 1 19691) . The capture of stars 
has been put forward as an explanation for ultra-luminous 



* We used their equation (3) rather than (4) and CLOUDY version 
05.07 rather than 07.02. 



X-ray sources (see e.g. lFabbianoll2006l . for a review). How- 
ever, we neglect this process in the current work, and instead 
focus on how BHs can accrete gas from their surroundings. 

The model presented in this section is a substantially 
modified version of the model introduced by SOS and em- 
ployed in almost all of the large-scale numerical simulations 
of AGN growth thus far available in the literature (see Table 
[2] for an overview). 

Because cosmological simulations have neither the res- 
olution nor the necessary physics to simulate the formation 
of the seed BHs that eventually grow into SMBHs, it is as- 
sumed that low-mass seed BHs are produced sufficiently reg- 
ularly that every halo above a certain threshold mass con- 
tains o ne such object at i ts cen ter. Here, our model follows 
that of iDi Matteo et "ahl l|2008l ) in that we regularly run a 
friends-of-friends group finder with linking length e qual to 
0.2 ti mes the initial mean inter-particle spacing (Davi s et al.l 
Il985h on all of the dark matter particles during the simula- 
tion. We do so at times spaced evenly in log expansion fac- 
tor, a, such that Aa = 0.02a, which corresponds to a proper 
time of ~250 Myr ('^ 70 Myr) at redshift zero (three) for our 
cosmology. When a halo grows above some threshold mass, 
JTihaio.min, and does not already contain a BH, then its most 
gravitationally bound baryonic particle is converted into a 
coUisionless BH particle. The initial mass of these BHs is 
usually chosen to be well below the resolution limit of our 
cosmological simulations (see Sec. |4]), and as such we need 
to employ sub-grid models to follow the BH. Although we 
convert the entire particle into a 'black hole particle', the 
mass of the seed BH (rrisccd) associated with this particle 
is usually initially significantly less than the particle mass 
(mseed rn^; where mg is the simulation gas particle mass). 
We therefore store the mass of the subgrid BH separately. 
For the gravitational interactions, other than BH accretion, 
the full mass of the particle (mg) is used, but for calculat- 
ing the BH-specific processes we use the sub-grid BH mass 
(?tibh). We now discuss in more detail the manner in which 
we track the growth of the BH. 

BH particles are coUisionless sink particles that con- 
tain a sub-grid BH, initially of mass mgecd, chosen to be 
well below the observed mass of BHs in haloes of this size. 
From their initial seed mass, BHs may grow via one of two 
processes: mergers with other BHs and accretion of sur- 
rounding ambient gas. We now treat each of these pro- 
cesses in turn. In our models BHs accrete from the sur- 
rounding ambient gas phase at a rat e proportional to that 
given by the Bondi-Hoy le-Lyttleton (|Bondi fc HovldHoij : 
iHovle fc Lvttletonlll93dh formula 



47vG menP 
/ 

■(ci +^2)3/2 



(1) 



where rriBH is the mass of the BH, Cs and p are the sound 
speed and gas density of the local medium, v is the veloc- 
ity of the BH relative to the ambient medium, and a is a 
dimensionless efficiency paramet er. The factor a did no t ap- 
pear in the origina l anal yses of iBondi fc HovQ l| 19441 ) and 
iHovle fc LvttletonI (|l939l ). but was introduced by SOS as a 
numerical correction factor, to compensate for the limita- 
tions of the numerical simulations. The assumption that BHs 
grow via Bondi-Hoyle accretion is reasonable even if they are 
in reality fed by accretion discs that are far smaller than the 
resolution limit of our simulations as long as the latter grow 
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by Bondi-Hoyle accretion. However, we will see that very 
large factors of ct are required for low-mass BHs to grow, in 
which case one cannot claim to be simulating Bondi-Hoyle 
accretion. 

The amount of accreted mass is related to the rate 
of growth of the BH bjjf| ttibh ~ fnaccr(l — tr), where 
is the radiative efficiency of a BH, which we always as- 
sume to be 10%, the mean v alue for the radiatively efficient 
IShakura k. SvunvaevI (119731 ) accretion onto a Schwarzschild 
BH. 

In order to resolve Bondi-Hoyle accretion onto a BH we 
need to reso lve the Bondi-Hoyle radius (rb), defined as (e.g. 
lEdgaiibOO^ I: 



rh = 



GiriE 



0.042 



(i^)(io^)"'^P" 



Comparing this to the Jeans length. 




(3) 



where Mj is the Jeans mass, we see that rb ~ Lj if rriBH ~ 
Mj and rb ^ Lj if ttibh S> Mj. Hence, any simulation 
that resolves the Jeans scales will also resolve accretion onto 
black holes of mass meH > rrig. We can then parametrize 
accretion onto a BH in two different ways: 

1. Density-Independent Accretion Efficiency: 
Most AGN models in the lite rature use a con- 
stant value of Q = 10^ (e.g. ISpringel et all l2005l: 
Pi Matteo et al.l l2005l : ISiiacki et all l2007l : iDi Matteo et"aLl 



2008': 'Bhattacharva et al. '2008': 'Colberg fc di Matted l2008l : 
Croft c t al..,200& : .Johansson et al. 2009.. see also Table [2L 
Although most authors do not motivate or even mention □ 
their choice of a, we note that values much greater than 
unity can be justified in one of two ways: firstly by not- 
ing that the Bondi-Hoyle accretion rate depends strongly 
upon the local ISM sound speed. Galaxy formation simula- 
tions currently have neither the resolution nor the physics to 
self-consistently track the properties of the cold- phase of the 
ISM, and as such the temperature of the gas accreted by the 
AGN may be overestimated by orders of magnitude. Hence, 
we can justify very large values of a in star forming gas. 
Secondly, in low-resolution simulations we do not resolve 
the Jeans scale, even in single-phase gas, so the density of 
the gas at the Bondi radius is underestimated, allowing us 
to again justify large values of a. We call models that use 
a fixed value of a 'constant-a' models. Constant-a models 
are parametrized by a constant multiplicative factor in the 



° Note that 305 neglected the (1 — e^) term and used mBH = 

^a ccr . 

^ ISprinpI et al.l ||2005|): iDi Matteo et al.l Jioosl'): ISiiacki et al.l 
ll2007l'): (pi Ma tteo ct a l.| bOO^ : iBhattacharva et al.l l|2008l) : 
IColbere fc di Matteg ^200^)all do not discuss or mention the 
value of a that they used, but we have been inforraed by 
V. Springel that they assumed a = 100. iKhalatvan et al~ 1 ll2008i ) 
state explicitly that in their models a = 300 and justify this by 
noting that when the density of the ISM is smoothed on the scale 
of the computational resolution, the recovered densities are much 
lower than would be ex pected on the sca l e of t he Bondi radius. 
Using similar reasoning, [Johansson et al ] jiooi) reach a similar 
conclusion and set a = 100. 



Bondi-Hoyle accretion rate, oq. An alternative way of in- 
creasing accretion rates is to employ a subgrid model for 
the unresolved ISM properties, and use this to artificially 
boost the ISM densities local to the BHs. We discuss this 
further in the following sections. 

We will show in Sec. 14.11 that the use of a constant-a 
model has a profound effect on the ability of BHs to grow, 
and that changes in this very poorly constrained parame- 
ter can lead to large changes in the global properties of the 
simulation such as the global density in BHs. We empha- 
size that values of a S> 1 imply the assumption that the 
simulation predictions for the gas density and temperature 
are sufficiently wrong that the Bondi-Hoyle accretion rate 
is underestimated by two orders of magnitude. Although 
this assumption can be justified for high-density gas, it does 
mean that the value of a is in fact more important than 
the predicted densities and temperatures. One could there- 
fore argue that models of this kind do not really simulate 
Bondi-Hoyle accretion. 

Cosmological simulations can, however, already model 
Bondi-Hoyle accretion of low-density gas and it hence makes 
sense to use an accretion model for which the "fudge factor" 
a becomes unity in the regime where the simulations are 
reliable. We therefore introduce a new class of black hole 
accretion models in which the the value of a depends on the 
local gas density, while keeping the number of free parame- 
ters fixed. 

2. Density-Dependent Accretion Efficiency: The 

assumptions used to justify large values of a in simulations 
similar to ours break down when two conditions are satisfied: 
firstly the local gas density must be lower than required 
for the formation of a cold (i.e. T <^ 10'' K) phase, and 
secondly the simulation must resolve the Jeans scale of the 
single-phase gas. Our highest resolution simulations (as well 
as many published AGN simulations) do resolve the Jeans 
scale at the star formation threshold and so in contrast to 
most published AGN models we choose to parametrize the 
accretion efficiency parameter as a function of density 



if riH < rig 
otherwise. 



(4) 



Here, the accretion efficiency (a) becomes unity for den- 
sities lower than the critical value required for the forma- 



tion of a cold interstellar gas phase 



= 0.1 cm- 



Sec. [2]). As discussed above, provided the simulations re- 
solve the Jeans scale, the Bondi radius will be resolved for 
BHs with JTiBH ^ rUg, which means that values of a ^ 1 
are unphysical for such BHs. We then choose to parametrize 
our lack of knowledge about both the physical properties of 
the multiphase ISM and the rate at which it accretes onto 
the central AGN using a power-law of the gas density, with 
slope /3. This constant-/? model, which has the same number 
of free parameters (one) as the constant-a models used in 
previous work, allows us to correctly describe accretion in 
the physical regime where it is resolved by our simulations 
and to introduce a reasonable scaling when it is not. We 
will show in Sec. I5.2l that the change from a constant-a to a 
constant-/? model can have a profound effect on the growth 
of BHs, particularly for low mass galaxies. We call models 
of this type 'constant-/? models'. 
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A second approach to boosting accretion rates, which 
operates in a similar manner to the constant-/? models, is 
to make use of a subgrid model for the unresolved sub- 
grid p hysics not encapsulated by the simulations. For ex- 
ample, iPelupessv eFaLl (|2007| ) use the star formation and 
supernova feedback models of SH03 to estimate the amount 
of time that a BH sp ends in dense, molecular clouds, and 
lOkamoto et"al] l|2008l ) use a subgrid model in which drag 
due to stellar radiation on a clumpy ISM can give rise to 
large accretion rates onto a central BH. We note that differ- 
ing implementations of the subgrid model can lead to large 
differences in the properties of the ISM and for the purposes 
of this work we emphasize that the functional form (Eq. |4]), 
as well as the value for P, are ad-hoc. Any function for which 
a — > 1 at gas densities for which the simulations are reliable 
and for which a 1 at higher densities would do. We chose 
to use a simple power-law dependence because it satisfies 
these constraints, is continuous and uses only one free pa- 
rameter. We will investigate the effect of changing /3 in the 
following sections. In the limit that /3 — > the behaviour of 
the constant-/? model will tend towards that of a constant- 
a model with qq = 1, and in the limit that /3 ^ oo the 
model tends towards behaviour where the accretion is pure 
Bondi-Hoyle in non star- forming gas, and always Edding- 
ton limited in gas with densities above the star formation 
threshold. We caution that this prescription is not suitable 
for simulations that resolve the relevant physics at densities 
exceeding njj. 

Because we have changed the density-dependence of the 
accretion rate, we cannot claim to be simulating Bondi- 
Hoyle accretion. Values of a 2> 1 are, however, motivated 
by the Bondi-Hoyle formula. Moreover, for tih > the 
density should be interpreted as the mass-weighted mean 
density of the unresolved, multiphase medium, smoothed on 
the scale of the spatial resolution of the simulation, whereas 
the density appearing in the Bondi-Hoyle formula applies to 
a single gas phase. Since the accretion rate-weighted mean 
density (which we can only compute if we know the mass 
distribution of the multiphase g function of density 

and temperature) is unlikely to be proportional to this ef- 
fective density, there is no reason to keep the Bondi-Hoyle 
scaling. For this reason, and because a S> 1 implies the as- 
sumption that the predicted densities and temperatures are 
greatly in error, we argue that constant-a prescriptions with 
Q S> 1 can no more claim to be modeling Bondi-Hoyle ac- 
cretion than constant-/? prescriptions. We prefer the latter 
since it allows us to get the right answer in the regime where 
the simulations are reliable, that is, at sufficiently low den- 
sities. Finally, we note that both accretion models use the 
Bondi-Hoyle scaling of the accretion rate with the mass of 
the BH, rhaccr oc m|jj. 

In common with the models of S05 we limit the accre- 
tion rate to the Eddington rate: 

ATvGmBiimp 
JTlEdd = , (5) 

where mp is the proton mass and ctt is the Thomson cross 
section for scattering of free electrons. Because rriBdd oc itibh 
whereas iriaccr oc jtibh i Eddington limited accretion tends to 
be more important for more massive BHs. 

Following SOS, we allow BH particles to stochastically 



swallow neighbouring baryonic particles with a probability 

_ f (mBH - mpart)p~^Vy(r-BH - Ti, /iBH ) if mflH > mpart 

otherwise 

where p is the local gas density, tubh is the mass of the sub- 
grid black hole, mpart is the mass of the particle containing 
the sub-grid BH and VF(rBH — ri, Hbh) is the SPH kernel, 
evaluated between the positions of the BH and gas parti- 
cle i. The BH smoothing length, /ibh, is chosen such that 
within a distance /ibh from the BH there are A^'ngb ~ 48 
neighbours, the same number of neighbours as we used in 
our SPH calculations. This process ensures that the mass of 
the BH particle always closely tracks ttibh 

When the mass of the BH particle is smaller than or of 
the same order of magnitude as the simulation mass resolu- 
tion, the black hole does not dominate the local dynamics 
and may wander from the centre of mass of its parent halo 
due to numerical effects. Conservation of momentum from 
accreted ISM gas can lead to similar effects. In order to avoid 
this we employ the same scheme as in the models of S05. At 
every timestep the gravitational potential energy is calcu- 
lated at the position of each of the BH's neighbouring gas 
particles and the BH particle is repositioned on top of the 
particle with the minimum potential energy. In order to pre- 
vent the BH from being 'dragged' by a minimum-potential 
particle with a large relative velocity, we only perform this 
process if the relative velocity between the BH and its most- 
bound gas particle neighbour is less than 0.25 Cs, where Cs 
is the local sound speed. This process ensures that the loca- 
tion of the BH particle always tracks the centre of mass of 
its parent halo very closely. This procedure is halted after 
the mass of the SMBH becomes greater than ten times the 
initial gas particle mass in the simulation because by this 
point the BH dominates the dynamics in the centre of the 
halo. 

3.2 Black hole mergers 

Galaxy mergers are thought to be one of the major processes 
driving the evolution of galaxies. When galaxies merge it is 
expected that their central BHs will eventually also merge. 
Indeed, the build-up of BHs through mergers may play an 
important part in the growth of SMBHs. Similarly to SOS 
we have implemented BH merging as follows. When any two 
BHs pass within a distance /ibh of each other with a relative 
velocity smaller than the circular velocity at a distance /ibh 
{vici < \/ GmBn/^BH, where /ibh and niBu are the smooth- 
ing length and mass of the most massive BH in the pair 
respectively) then they are allowed to merge. This velocity 
criterion is necessary in order to prevent BHs from merging 
during a fiy-through encounter of two galaxies, as this could 
lead to BHs being quickly removed from their host galaxies 
due to momentum conservation. This velocity scale is some- 
what different from that employed by SOS, who used the 
local sound speed, Cs, as the relevant velocity scale, argu- 
ing that the sound speed represents a simple measure of the 
characteristic velocity scale of the galaxies, and hence gives 
a simple measure of the velocity scale at which BHs will be 
able to merge. However, because AGN input large amounts 
of energy into their surroundings, it is not necessarily true 
that the sound speed local to the AGN reflects the depth of 
the potential well. 



Simulating SMBH growth and feedback 7 



The BH merging rate estimated from our simulations 
likely represents an upper limit to the true merger rate as our 
simulations do not have the resolution required to resolve the 
formation of the tight BH bin aries that are a prer equisite 
for their eventual coalescence (ICallegari et aljboosl ). Since 
it is not yet fully understood how lo ng it takes to harden 
a BH binary l|Makino fc Funatdliool) . we assume that the 
merging process is instantaneous. 



3.3 Energy feedback from black holes 

The precise mechanism by which energy emitted from a 
BH is coupled to the surrounding medium is as yet un- 
known, but plausible mechanisms include radiation pres- 
sure on free electrons (which gives rise to the classical Ed- 
dington limit), Compton heating of the infalling gas (e.g. 
ICiot ti fc Ostrikcr 2001; Wang et al. 2006), photoionization 
pressure (Buff fc McCrav ,1974 : .Cowie et al. 1978) and ra - 
diation pressure on dust grains fe.g. iMurrav et all |2005| ). 
Regardless of the precise coupling mechanism, there is a cat- 
alogue of observational evidence indicating that energy out- 
put from AGN can drive galactic outflows. For e xample, ab- 
sorbe rs seen in X-rays show evidence of outflow (jLaor et al.l 
1 19971 ) and broad absorption line s ystems show eviden ce of 
outflows at very high velocity (e.g. [Pounds et al.ll2003h . Al- 
though these observations indicate that high velocity out- 
flows are present around some AGN, they do not tell us how 
much mass (and hence how much energy) is present in the 
outflow. Estimates of the mass outfl ow rate in the wi nds are 
highly uncertain. Some studies (e.g. ICheloucIi3l2008l ) imply 
that the actual rate of mass outflow is only a small frac- 
tion of the bolom etric luminosity of the AGN sources, while 
other studies (e.g. lArav et al ]|2002l .[ 20081 ) suggest large mass 
outflow rates in quasar driven winds. 

In our models BHs inject a flxed fraction of the rest 
mass energy of the gas they accrete into the surrounding 
medium. The feedback is implemented thermally, that is: 
energy is deposited into the surrounding gas by increasing 
its internal energy, as opposed to the kinetic feedback used 
to inject supernova energy, which is deposited by kicking 
the gas particles (see Sec. [2} . The fraction of the accreted 
rest mass energy that is injected is assumed to be indepen- 
dent of both the environment and the accretion rate. We 
thus do not differentiate between 'qu asar mode' and 'ra - 
dio mode' feedback as in the models of lSiiacki et al.l (|2007h . 
In a future work we will consider how spatially distributed 
AGN heating mechanisms affect the cosmological evolution 
of galaxies. The amount of energy returned by a BH to its 
surrounding medium in a timestep At is given by 

-Efccd = CfermBHC^At , (6) 

where ef is the efficiency with which a BH couples the radi- 
ated energy into its surroundings - a free parameter in our 
simulations - and c is the speed of light. Only the product 
of Er and Ef is important in calculating the amount of energy 
feedback in our model. 

Because our sub-grid model for SF relies on an effective 
EOS and does not include a (semi-)analytic sub-grid model 
for the multiphase ISM, our energy distribution mechanism 
is different from that in SOS. In contrast, because we pre- 
fer to minimize the use of semi-analytic models within our 



hydrodynamical simulations, our models rely only on an ef- 
fective EOS and leave the distribution of the mass over unre- 
solved gas phases undefined. We therefore need to make two 
changes to the EOS model of Schave fc Dalla Vecc hia (20C)J) 
that was used in our simulations without AGN feedback. 

Firstly, in the original models, once gas was identified 
as star- forming it was forced to remain on the EOS, un- 
til its density dropped below the critical density for star- 
formation, tih, it turned into a star particle, or it was kicked 
into the wind. It is therefore necessary that we change this 
by allowing strongly heated gas to leave the EOS. This is im- 
plemented numerically by taking gas that is heated by more 
than 0.5 dex above the EOS in a single time step off the EOS 
(i.e. it is no longer star-forming and its pressure is no longer 
constrained to lie on the EOS). Gas is placed back onto the 
EOS if its temperature falls back below 0.5 dex above the 
EOS temperature corresponding to its density. By checking 
SFRs, both globally and for individual objects, and by com- 
paring gas distributions on the p — T plane we have verified 
that making this change to our EOS model has a negligible 
effect on the results in a simulation that does not include 
AGN feedback. A second possible change to the AGN model 
would have been to treat the effective EOS as a lower limit 
to the gas temperature. We tested this and again found the 
differences in our results to be negligible. We choose to use 
the first procedure in order to facilitate direct comparisons 
between the simulations containing AGN and those run ear- 
lier in the project. 

Secondly, in order to ensure that the thermal feedback 
from BHs is not immediately radiated away it is necessary 
to impose a minimum heating temperature. BHs store feed- 
back energy until they have accumulated enough energy to 
increase the temperature of Whcat of their neighbours by an 
amount of ATmin, i.e. 

J-, _ nheatmgfcBArmin 



where -Edit is the critical energy for a heating event to be 
triggered and is the mean molecular weight of the gas 
(we assume /i = 0.58, appropriate for a fully ionized gas of 
primordial composition). The internal energy of the heated 
gas is instantaneously increased by an amount i?crit. This 
implementation of quasar mode fee dback is simi l ar to the 
radio mode feedback introduced by ISiiacki et all l|2007l ). If 
ATmin is set too low then the cooling time of the AGN 
heated gas will remain very short, and the energy will be effi- 
ciently radiated away. If riheat ATmin is set too high then the 
threshold energy for a heating event to occur and hence the 
time period between AGN heating events will become very 
large. In particular, a time interval larger than the Salpeter 
time would prevent the BH from regulating its growth. Fi- 
nally, we note that the energy is deposited into the ambient 
gas isotropically, equally distributed to a random fraction 
TT-iioat/A^ngb of the BH's neighbours. If, on a given timestep, 
a BH accretes more energy than necessary to heat nhcat par- 
ticles to ATmin then the process it repeated until the BH has 
distributed all of its energy, so individual gas particles may 
be heated by an amount ATnin multiple times on a given 
timestep. 
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4 PARAMETER CHOICES 

Both the mechanism by which BHs grow and the efficiency of 
their thermal feedback can be changed drastically by chang- 
ing the values of the parameters of the AGN model. In this 
section we discuss how parameter values are chosen to min- 
imize unphysical numerical effects whilst simultaneously re- 
quiring that the global properties of the BH distribution sat- 
isfy various observational constraints. For quick reference, 
Table [1] contains a full list of the parameters that control 
the behaviour of the BH growth and AGN feedback model, 
along with their fiducial values. 

Because it is difficult to discuss the effect of each pa- 
rameter in isolation, we will first discuss the general prop- 
erties of the BH model (e.g. growth mechanisms, feedback 
efficiency) and use each of these general themes to motivate 
our fiducial choices for the parameters of the AGN model. 




-3 -2 



-1 



log(n [cm"^]) 



4.1 Black hole growth 

We allow BHs to grow by two processes: mergers with other 
BHs and accretion of gas. The BH accretion time-scale 
taccr = 771bh/"1bh IS, for Boudi-Hoylc accretion (Eq. 
proportional to m^^. Therefore, depending upon choices for 
various parameters, the initial growth of BHs can proceed 
in one of two different ways. Firstly, in the regime where 
the time scales over which gas accretion operates are very 
long, BHs may grow primarily by mergers with other (seed 
mass) BHs until rriBH is large enough for the accretion rate 
to become appreciable. In this regime black holes initially 
grow at a rate governed by the integrated mass of seed BHs 
they collide with. Secondly, seed mass BHs may have accre- 
tion rates large enough for the BHs to experience runaway 
growth until their accretion rate is limited by feedback pro- 
cesses. 

We can estimate the growth rate of BHs by noting 
that in our star formation model we impose an effective 
equation of state on star-forming gas, and as such can im- 
mediately calculate the local ISM pressure, P, (and hence 
Cs — VtP/p) from 

p=p.4^r\ (8) 

where njj = 0.1 cm~^ and Pcrit are the critical threshold den- 
sity and pressure for star formation respectively (see Sec. [2]). 
Fig. [T] shows BH growth times as a function of the ambient 
gas density for both constant-a and constant-/? models, as- 
suming here that the BH is of mass 10^ M0. For the pur- 
poses of this plot we assume that when gas densities are 
below the star-formation threshold the EOS is isothermal, 
but note that in the simulations we calculate the pressure 
self-consistently. Following other authors, we set ao = 100 
for the constant-Q model. For the constant-/? lines we set 
P = [1,2,4]. The horizontal, red line in this plot shows 
the growth time of a BH that is accreting at the Edding- 
ton rate (i.e. the Salpeter time). The Salpeter time depends 
only upon physical constants and the BH radiative efficiency, 
such that 

isaipotor = -. = = 4.5 X 10 hr-r yr • (9) 

yriEdd ATvGmp VO.l/ 

It is immediately clear from Fig. [T]that the choice of accre- 
tion model strongly affects the local density at which the BH 



Figure 1. BH growth times (niBH/"^BH) a function of the 
ambient gas density under various accretion models, all for a BH 
mass of 10® Mq. The normalization of all black lines scales as 
rrigjlj. Lines are shown for both a constant-a accretion model 
(solid, black/grey line) and for constant-/? accretion models (all 
other black/grey lines). The solid, red line shows the Salpeter 
time (the growth time for a BH accreting at the Eddington rate), 
and represents the lower limit on the BH growth time in the 
simulations. The grey section of each line represents the region 
where the accretion rate is greater than the Eddington rate. The 
vertical dotted line shows the star formation density threshold, 
= 10~^ cm~^. Above this density the gas follows the effec- 
tive equation of state defined by Eq. \E\ For lower densities we 
have assumed the gas to be isothermal (only in this figure, not 
in the simulations). Note that the constant-Q accretion model 
predicts that a 10® Mq BH will be growing at an Eddington lim- 
ited rate even in gas with density at the star-formation threshold 
(0.1 cm-3). 

growth becomes Eddington limited, with black holes accret- 
ing in the constant-a model becoming Eddington limited at 
densities 1-2 orders of magnitude lower than the same black 
hole accreting in the constant-/? model. 

From the simulations we find that, for our chosen value 
of rrihaio.min (^IOOttidm; See Sec. 14. 3p and at high redshift, 
typical birth densities of BHs are ~ 10 — 100 times the star 
formation threshold. Hence, in the regimes of interest, in 
constant-a models all BHs of mass > 10^ M0 grow initially 
at close to the Eddington rat^ll (for a seed mass of 10^ Mo 
and typical initial gas densities of 10^ — 10^ cm~'^ the initial 
Eddington ratio is 0.037-0.37) until feedback effects reduce 
the local gas density to values below the threshold for star 
formation. From Fig. [1] we can see that for a BH of mass 
lO'' M0 the accretion rate remains Eddington limited until 
the local gas density falls to nn < lO"'^ cm~^. Fig.[2]presents 
this information in a slightly different manner by showing 
the density above which a BH's accretion rate is Eddington 
limited as a function of BH mass for the same accretion mod- 

^ Note that the choice of sub-grid ISM model is a significant 
source of additional uncertainty in the accretion rates. For exam- 
ple, the use of the SH03 sub-grid multi-phase ISM can give rise 
to differences of almost an order of magnitude in the accretion 
rates of seed mass black holes due to the use of different effective 
equations of state. 
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Table 1. Pull list of parameters of the AGN model, along with the values they take in the fiducial set of simulations, and short definitions. 
Parameter Piducial Value Description 



^sccd 

^halo,min 

fir 

ao 



0.001 rrig Initial mass of the sub-grid BH 

100 rriDM Minimum halo mass into which BH seeds may be placed 

0.1 Radiative efficiency of the BH accretion discs 

0.15 Praction of energy emitted by BHs that couples into the ambient gas 

1 Number of neighbouring particles heated per feedback event 
10* K Amount by which BH feedback heats surrounding gas 

Depending on the accretion model, one of the following parameters is used: 

100 Normalization of the Bondi-Hoyle accretion efficiency (Pq. [ij in constant-a models 

2 Slope of the Bondi-Hoyle accretion efficiency (Eq. |4]l in constant-f) models 




4 6 8 10 
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Figure 2. The gas density above which the accretion rate onto a 
black hole becomes Eddington limited as a function of BH mass, 
assuming that gas with density above the critical density for star 
formation, njj = 10~^ cm~^ has properties governed by the ef- 
fective EOS and that gas below the critical density follows an 
isothermal EOS. The thick, black line shows the behaviour of a 
constant-a accretion model, and all other lines show how models 
with a constant-/? accretion rate behave. The grey line shows the 
critical density for star formation in our simulations. Except for 
very low BH masses, the constant-a model becomes Eddington 
limited at much lower gas densities than the constant-/? models. 

els as in Fig. [1] Again, it is clear that the gas density below 
which the accretion rate depends on the density, and can 
thus more easily be regulated by feedback from the AGN, 
has a strong dependence on the accretion model used. Take, 
for example, the case of a BH of mass 10* M© , here in a 
constant-a model this BH will grow at the Eddington rate 
until it can modulate its local density to below 10~^'^ cm~"^, 
more than two orders of magnitude below the star-formation 
threshold. This is well within the regime where Bondi-Hoyle 
accretion is resolved in the simulations. 

In the constant-a model with ao = 100 the growth of 
BHs therefore proceeds as follows: Seed mass BHs (typical 
seed masses are in the range 10"^ — IO^Mq in our simula- 
tions) grow exponentially by Eddington limited accretion 
until feedback from the BH has decreased the local ISM 
density to the point that growth is no longer Eddington 
limited, and further energy output from the AGN can de- 
crease the accretion rate. For BHs with masses greater than 
10^ Mq self-regulation can only occur at densities orders of 



magnitude below the star formation threshold (Fig. [2)). In 
this regime we resolve Bondi-Hoyle accretion, invalidating 
the assumption used to justify large values of a in the first 
place. 

In contrast, simulations that employ a constant-/? 
model, are not necessarily Eddington limited from birth. 
Taking again the case of a 10^ BH, and our fiducial value 
of /? = 2, we see that the BH can decrease its accretion rate 
at a much higher gas density, and as such the period of Ed- 
dington limited growth will be much shorter than for the 
constant-a model. The difference in the gas density below 
which AGN accretion rates are no longer Eddington limited 
can lead to large differences in the properties of low mass 
galaxies and BH growth in small haloes. 



4.2 Efficient thermal feedbaclc 

As discussed in Sec. 13.31 it is necessary to impose a minimum 
heating temperature in order to prevent BHs from heating 
their surroundings before they have generated enough en- 
ergy for thermal feedback to become efficient. 

Two parameters control how efficient this feedback may 
be: the number of neighbours to heat (rihoat) and the tem- 
perature to which the neighbours are heated (ATmin). Two 
competing effects control our choices for these parameters. 
If ATmin is too low thcu AGN heated gas will retain a low 
temperature and therefore also a short cooling time (an anal- 
ogous problem to the 'overcooling' of supernova hea ted gas 
in early cosmological simulations (|Katz et al.lll996l ). which 
is, however, usually attributed to an overestimate of the gas 
density) . In this regime the energy will be immediately radi- 
ated away, making AGN feedback ineffective. Conversely, if 
ATmin or nhoat are set too high then the time scale over 
which BHs accrete enough energy to heat riheat of their 
neighbours by an amount ATmin will become longer than 
either the dynamical time in the vicinity of the BH (or the 
Salpeter time for Eddington limited growth), leading to spu- 
rious growth as BHs are unable to self-regulate. 

The choice of the minimum heating temperature, 
ATmin, is motivated by the fact that we wish to choose 
the minimum value (and so minimum time between heating 
events) for which the cooling time of the heated gas is long 
enough that the energy is not immediately radiated away. In 
practice it was found that ATmin = 10* K is the minimum 
temperature for which BH feedback has a sufficient effect on 
galaxy clusters. We return to this point in Sec. 15.21 

The second parameter, njioat, is calibrated by noting 
that although ideally we would like to allow AGN to heat 
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gas instantaneously, our finite resolution forces us to store 
energy until the feedback is effective, hence introducing a 
delay to AGN heating. We can minimise the effect of this 
delay by noting that the numerically imposed time between 
heating events should be lower than the typical time-scales 
of dynamical processes that affect AGN feedback. If nheat 
is set too high then it is possible that the amount of time 
taken for a BH to accrete enough energy to perform a heat- 
ing event would be large enough that we see spurious growth. 
We can quantify this effect by calculating the mean time be- 
tween heating events for BHs of different masses in different 
density environments. This is demonstrated in Fig. [S] which 
shows the mean time between heating events as a function 
of msH/jTig for models with both constant-/? and constant-a 
accretion rates. Plotting the heating time as a function of 
this mass ratio allows us to make this plot in a resolution 
independent manner. Fig. [3] assumes ATmin = 10* K and 
Jiheat ~ 1, but all lines can be shifted vertically in propor- 
tion with the quantity ATmin?iheat . To make the time be- 
tween feedback events as small as possible we should choose 
JT-hoat as small as possible, but it is not immediately obvious 
that BHs will be able to regulate their growth if we heat 
only a small number of neighbours of an AGN that it will 
be able to self-regulate. However, we found from numerical 
tests (Sec. 15. 2|) that nhcat = 1 is sufficient for BH feedback 
to be effective, and so riheat = 1 is the parameter value used 
in our fiducial simulations. The dependence of our results 
on the two parameters ATmin and rihcat will be discussed in 
Sec. [521 

The next model parameter is ef, the efficiency with 
which energy radiated from the BH is coupled to the ISM. 
The parameter ef sets the normalizations of the global BH 
density and the BH-galaxy scaling relations. We therefore 
tune ef after setting all of the other parameters in order to 
match the redshift zero observations (Sec. 15.21 Fig. [7jl) and 
find that a value of ef — 0.15 provides a good match to the 
observations. 

4.3 Black hole seed mass and minimum halo mass 

Our initial choice for the halo mass into which we insert a 
BH is motivated by the fact that we wish for every resolved 
halo with rrihaio ^ 'mseed to contain a seed BH. We therefore 
choose to place BH seeds into haloes of a constant particle 
number. Using mhaio,min = 100 muM ensu res that haloes 
containing BHs are always well defined (e.g. iDiemand et al.l 
I2OO7I ). The choice of a constant particle number halo mass 
also has the advantage that if we change the simulation mass 
resolution, BHs will still be placed into the smallest allow- 
able mass of dark matter haloes without the need to tune 
any parameters. Note, however, that this prescription will 
have to be changed for simulations that have sufficient res- 
olution for 100 moM to be comparable or smaller than the 
minimum halo mass expected to be able to form seed mass 
BHs. 

Given the minimum halo mass into which we place BH 
seeds, we must ensure that the integrated number of seed 
BHs generated between redshifts z — 00 and zero is much 
smaller than the observed cosmic BH density. We can obtain 
an upper limit on the cumulative cosmic density of BH seeds 
by taking the redshift zero dark matter halo mass function 
f{m) — n(m)dm assuming that all collapsed mass was as- 



sembled through mergers of critical mass haloes: 

/5socd(mhaio,min) < / mf{m)dm. (10) 

,IIlin J JTii-,-.],-, i-niri 

This quantity is plotted for a number of values of the seed 
BH mass in Fig.|4l where the two vertical grey lines represent 
the masses of haloes of 100 DM particles (=mhaio,min) in our 
fiducial simulations run at the mass resolution as the OWLS 
runs of Schaye et al. (in preparation). 

Given choices for ATmin and nhcat, we can use Fig. [3] to 
place a minimum limit on the BH seed mass, rrisccd- Here, we 
show the time between heating events as a function of BH 
mass, for both constant-a and constant-/? models. The time 
between heating events for BHs accreting at the Eddington 
accretion rate are shown as red lines in each panel. We now 
note that we require BH heating to occur regularly in high 
density environments. In particular, in order for a BH to be 
able to effectively self-regulate its own growth, we require 
that the numerically imposed minimum duty-cycle, Athoat, 
is less than the Salpeter time, the characteristic growth time 
for black holes accreting at the Eddington rate. It is clear 
from Fig.Oby comparing the BH Salpeter time (blue line) to 
the BH duty cycle that in high density environments (nn > 
10^ cm""^) this condition is satisfied only if the BH mass is 
greater than 10~^mg. This provides a minimum allowed seed 
mass in our models. 

However, in addition to ensuring that BHs grow in a 
physical manner and that their feedback can be effective, 
we must also satisfy various observational constraints. Most 
fundamentally, it is known that the present day cosmic BH 
density is (4.2±1.1) x lO'^ Mp /Mpc^ jsh ankar et "al]|2004l ') or 
{4:.22till) X lO'' Mo/Mpc^ (|Marconie t al. 2004j), although 
we caution that a more accurate consideration of the ef- 
fects of cosmology may lead to a slightly higher d etermina- 
tion of the BH density (|Graham fc Driveill2007bl ). In order 
that we do not violate these observational constraints in the 
presence of substantial BH growth through accretion, we re- 
quire the z — global seed density to be much smaller than 
the observed BH mass density. We now ensure that - given 
all of our other parameter choices - rrigcod = 10~^ rrig does 
not violate this constraint on the global BH density. For 
our simulations 10"'^ rUg corresponds to a BH seed mass of 
1.2 X 10^ Mq. It is clear from Fig.|3]that the maximum pos- 
sible contribution of the seed BH mass to the cosmic density 
is at least a factor of 10 less than the redshift zero obser- 
vations, which we indicate by the grey, horizontal shaded 
region. We will see in Fig. [7] that, as expected, the actual 
contributions of seed BHs to the total cosmic density are 
much smaller than this value. 

An additional observational constraint is placed by the 
well-defined relation between the mass of a BH and the mass 
of the bulge compone nt of a galaxy, uibh ~ 0.006 mbuigo 
(|Magorrian et alll 19981 ) . In simulations it is more convenient 
to work with the relation bet ween BH mass an d dark mat- 
ter halo mass, investigated bv iFerraresd (|2002l ). who found 
that in halos of 10^^ Mq the ratio mBH/niha.io ~ 10~^. This 
provides a second, related constraint on the mass of seed 
BHs: we wish to place them below this relationship so that 
they can subsequently grow on to the observed redshift zero 
relatiorlfl For our parameter choices (mseed = 10~^ nig and 

* In some BH seed generation scenarios, for example the direct 
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Figure 3. Time between AGN heating events for Bondi-Hoyle (black and grey lines) and Eddington-limited (red lines) accretion as a 
function of the ratio between the BH mass and the simulation gas particle mass. The grey section of each line represents the region 
where m > rh^^^. In this plot we assume that riheat = 1 and ^T^i^ = 10* K. All lines may be shifted vertically in proportion with 
nhcat ATlnin . The horizontal blue line shows the Salpeter time (mBu/'TtEdd)- Left panel: A constant-o accretion model with ao = 100. 
Right panel: A constant-/? accretion model with (3 = 2. The spacings between the lines change when nn > n^*, because the accretion 
efficiency becomes proportional to nj^. For any Eddington limited accretion and BH masses greater than 10 rrig the time between 
heating events is shorter than the Salpeter time, enabling the BHs to self-regulate. 



JTlhalo.min = 100 mDM) 
"Isccd 10~^ 



^halo,min 100 Vn^ - !^ 



■ 2.1 X 10" 



(11) 



where the last equality assumes our chosen cosmology. This 
ratio is indeed much smaller than the observed value. 

Finally, we note that our fiducial choice of seed mass, 
niseed = 10~^mg, will need to be modified for simulations 
that have sufficient resolution for this value to be below the 
expected BH seed masses. 

4.4 Comparison with previous work 

Through the arguments in the previous sections we were 
able to specify values for all of our model parameters. These 
fiducial parameter values are summarised in Table [1] The 
AGN model developed in this paper is a modification of that 
introduced in SOS, and used thereafter in a large number of 
works. As such it is instructive to compare our parameter 
choices with those employed in other studies, as collected in 
Tabled 

We turn our attention first to the AGN feedback effi- 
ciency, Ef. The value used in the present study (ef = 0.15) 
is significantly higher than that used in previous published 
studies, which all assume ef = 0.05 — 0.1. We can account 
for this difference if we note that, unlike the other studies, 
we do not employ the SH03 subgrid model for the ISM. Use 



collapse of matter in haloes we expect BH seeds to reside initially 
above the observed scaling relations. We will show in Sec. 15.2.21 
that in our models BHs grow onto the BH scaling relations re- 
gardless of whether they are initially placed above or below the 
relations. 



of a different subgrid model for the unresolved ISM is likely 
to lead to differences in the amount of radiative losses, as 
the effective density and temperature of the ISM differ sig- 
nificantly between the two models. 

We note that apart from differences in the ISM model 
and AGN heating mechanisms, the strength of the AGN 
feedback depends only on the parameter combination eftr, 
and that there is significant leeway in the value of e^. All 
studies presented in Table [2] assume — 10%, but val- 
ues close to tr = 20% are p ossible for thin-disc accr etion 
on to a Kerr BH (|Yu fc Trem ainc 20o3: lThorii3ll974l ). Re- 
cent observational determinations of span the fuU rang e 
of allowabl e values: = 3 - 35% JWang et al.1 |2006| ). 
e = 15% (lElvis et al.ll2002l: IYu fc Lull2008l), er = ^7 - 8% 
(|Cao fc Lill2008l : iMartinez-Sansigre fc Tavlodfioosl ) depend- 
ing upon the specific assumptions and models used in each 
study. 

The ratio of the minimum halo mass to the seed mass 
is similar in all of t he cosmological studies , with the excep- 
tion of the work of ^Khalatv an et al] (|2008^ , who performed 
zoomed cosmological simulations of an individual object. In 
order to avoid numerical issues when mBH < rrig, these au- 
thors forced the BH to accrete very quickly at early times 
before artificially halting its accretion until the stellar mass 
of the halo becomes large enough that the BH lies on the 
observed men — rn^, relation. 

We see also from Table [2] that the minimum halo mass 
(Tthaio,!!!!!!) is consistent between all of the published cos- 
mological studies to within a factor of 5 (1 — 5 x 10^" Mq). 
We will show in Sec. 15.21 that, for the constant-a accretion 
model assumed in previous studies, AGN feedback affects 
all haloes into which seed mass black holes are placed, and 
that changing mhaio.min by a factor of of ten has a large ef- 
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Figure 4. Maximum possible contribution to the cosmic BH den- 
sity from seed mass BHs, as a function of tlie minimum dark mat- 
ter halo mass, mjiaio mini assuming the dark matter mass function 
of Reed et al. (2006). Each black line corresponds to a different 
BH seed mass, as indicated in the legend. The horizontal, grey 
shaded region shows the observed cosmic black hole density at 
redshift 2 = (Shankar et al. 2004), and the vertical line indi- 
cates mhalo min for our fiducial simulation, which has a DM mass 
resolution of 8.64 X 1 — Mq. For our fiducial BH seed mass of 
rriseed = 1-2 X 1 — ^ M0 we see that the maximum possible contri- 
bution of seed mass BHs to the global BH density is much lower 
than the 2 = observations. 



feet on the BH properties of the less massive haloes in the 
simulation. Care should therefore be taken when comparing 
results of different simulations. 

The two areas where our work differs most from previ- 
ous models are the Bff accretion model and the SN feedback 
model. We turn our attention first to the SN model and note 
that almost all previous studies e mploy the work of SIf03, 
whereas we use the method of iDalla Vecchia fc Schavel 
(2008'). Contrary to SH03, SN winds in our model are local 
and not hydrodynamically decoupled from the surrounding 
gas. Hydrodynamical decoupling of supernova heated gas, as 
implemented in the SH03 model, guarantees that when gas 
is kicked it is able to escape the ISM (although it may sub- 
sequently return). On the other hand, if the hydrodynamic 
forces are taken into account, the gas can remain confined 
by pressure forces, and if it does manage to escape it may 
drag along much of its neighb ouring gas. It was shown in 
IPalla Vecchia fc Schavel l)2008h that in high-resolution sim- 
ulations of individual disc galaxies this change fundamen- 
tally alters the structure of the galactic disc and that hydro- 
dynamically coupled winds and generate galactic outfiows 
with properties broadly comparable with observations. We 
demonstrate in a companion work that SN feedback has a 
large effect on the properties of the AGN population and, as 
such, it is important to investigate a number of SN feedback 
prescriptions. 

The second area where our model differs significantly 
from others in the literature is in the accretion model. In the 
nomenclature of this paper, all of the previous AGN models 
are constant-a accretion models in which ao = 100 — 300. We 
show in later sections that the accretion model represents 
one of the most crucial elements of an AGN model, and 



that all results are very sensitive to the way in which BHs 
are allowed to accrete. 

We show in this work that the results derived from 
simulations depend on aspects of the BH model that are 
homogeneous between the studies that have thus far been 
published. The present work - which is carried out using 
different techniques and parametrizations for much of the 
sub-grid modelling - therefore provides a way to investigate 
the robustness of the models. 



5 SIMULATION RESULTS 

In this section we first introduce the simulation set used in 
this paper (Sec. 15. ip before demonstrating that the fiducial 
simulations reproduce redshift zero observational results and 
quantifying how robust our model is to poorly constrained 
parameter choices (Sec. 15. 2p . To begin with, however, we 
show for illustrative purposes the gas density in a 3 Mpc/fe 
thick slice from our 100 Mpc//i (LIOOA'^512) simulation at 
redshift zero (Fig. [SJ. Each panel represents a successive 
factor of five zoom. The bottom-left panel shows a region 
800 kpc/h across, centered on a 3 x 10^ M0 BH, contained 
in a halo with a stellar mass of 3 x 10^" M0 and a dark 
matter mass of 2 x lO^'^ Mq . Circles in this plot represent 
the locations of BHs with the area of each circle proportional 
to the logarithm of the mass of the BH. 

5.1 Simulation list 

In order to explore parameter space we have run a large 
number of smaller simulations, the details of which are 
summarised in Table [3] Simulation names are of the form 
LxxxNyyy, where xxx represents the simulation box size in 
comoving Mpc//i and yyy is the cube root of the initial num- 
ber of dark matter and gas particles. For example, the sim- 
ulation denoted I/100A''512 refers to a comoving simulation 
volume of 100 Mpc/h, which contains 512^ dark matter par- 
ticles and an equal number of baryonic particles. Simulations 
for which one of the parameters was changed from its de- 
fault value are denoted by appending a descriptive suffix to 
the end of the simulation name. For example, simulations 
without AGN feedback are named LxxxNyyy NOAGN, and 
correspond to the simulations denoted REF_LxxxNyyy in 
the OWLS project (Schaye et al., in preparation). 

5.2 The effects of changing the model parameters 

We now turn our attention to the effect of varying each of 
the parameters of our AGN model away from those selected 
in Sec. jj] In order to quantify the effect of different aspects 
of the AGN model on the results of our simulation, we split 
the model parameters into four separate categories: the ac- 
cretion model (ao; 0); the seed generation model (mhaio,min; 
JTisoed); the feedback efficiency (ef); and the heat distribu- 
tion model (ATmin; Jihoat). We look separately at the effects 
of changes in each of these parameter sets, and additionally 
consider two purely numerical effects: the simulation mass 
resolution and box size. For each set of simulations we make 
four diagnostic plots: in Fig. [6] we show the cosmic SFR 
density as a function of redshift; Fig. [7] shows the evolution 
of the global BH density, and the cumulative BH density 
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Table 2. Model parameters for AGN feedback studies in the literature. Columns are as follows: 1) Reference; 2) Efficiency with which 
energy emitted by a BH is coupled into the ambient gas; 3) Radiative efficiency of BH accretion discs; 4) Multiplication factor for the 
Bondi-Hoyle accretion rate; 5)For star forming gas the Bondi-Hoyle accretion rate is multiplied by (nu/lO"^ cm"'^)'^ (Eq. IIJ; 6) Mass 
of seed BHs; 7) Minimum halo mass in which black hole seeds are placed; 8) Number of DM particles corresponding to a halo of mass 
mjiaio mini ranges indicate the difference between the highest and lowest resolution simulations used in each study; 9) The ratio of the BH 
seed mass to the minimum halo mass; 10) Type of simulation: Iso., isolated model disc galaxy. Zoom, zoomed simulation of individual 
object in a cosmological volume. Cosmo., uniform mass resolution simulation of a cosmological volume; 11) Reference for star formation 
model used in the study: Springel & Hernquist (2003) (SH03), Schaye & Dalla Vecchia (2008) (SD08); 12) Reference for supernova wind 
model used in this study: Springel & Hernquist (2003) (SH03), Dalla Vecchia & Schaye (2008) (DS08) 



Study 
(1) 


(2) 


(3) 


«0 

(4) 


/3 
(5) 


™socd 
Mq 

(6) 


™halo,min 
Mq 

(7) 


-^halo,miii 

(8) 


"Isood 
"Ihalo.min 

(9) 


Type 
(10) 


SF Model 
(11) 


Wind Model 
(12) 


Springel et al. (2005) 


0.05 


0.1 


100 





10^ 


n/a 


n/a 


n/a 


Iso. 


SH03 


SH03 


Robertson et al. (2006) 


0.05 


0.1 


100 





10^ 


n/a 


n/a 


n/a 


Iso. 


SH03 


SH03 


Sijacki et al. (2007) 


0.05 


0.1 


100 





10^ 


5 X lO^O 


2260 


2 X 10"*^ 


Zoom 


SH03 


SH03 


Sijacki et al. (2007) 


0.05 


0.1 


100 





10^ 


5 X lO^O 


285-606 


2 X 10"'^ 


Cosmo. 


SH03 


SH03 


Johansson et al. (2008) 


0.05 


0.1 


100 





10^ 


n/a 


n/a 


n/a 


Iso. 


SH03 


SH03 


Di Matteo et al. (2008) 


0.05 


0.1 


100 





10^ 


1 X lO^O 


36-363 


1 X 10-5 


Cosmo. 


SH03 


SH03 


Khalatyan et al. (2008) 


0.1 


0.1 


300 





n/a 


n/a 


1463 


n/a 


Zoom 


SH03 


SH03 


This study 


0.15 


0.1 


1 


2 




lOOmDM 


100 


2 X 10"^ 


Cosmo. 


SD08 


DS08 



Table 3. Summary of simulation parameters. Columns are as follows: 1) Name of simulation; 2) Comoving box size; 3) Number of 
DM and gas particles; 4) Final redshift; 5) Multiplication factor for the Bondi-Hoyle accretion rate (Eq. [TJ; 6) For star forming gas the 
Bondi-Hoyle accretion rate is multiplied by {nii/10~^ cm~^)l^ (Eq. IIJ; 7) Fraction of radiated energy that is coupled back into the ISM 
(Eq. [SJ; 8) Initial mass of sub-grid BH; 9) Minimum halo mass into which a seed BH is placed; 10) Minimum AGN heating temperature 
difference; 11) Number of particles heated per event. Entries in boldface represent parameters that have been changed from the value 
they take in the default model. 

Identifier ^ iV.as ^end «0 /3 ^ ^^^^^^ n^..t 

(1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) 
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0.001 
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1 
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,0 
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1 


2 
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0.001 


100 


1 
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,0 


256^ 





1 


2 
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0.001 


100 


1 


1 


L050N256VHIEPS 
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,0 


256^ 





1 


2 


0.6 


0.001 


100 


1 


1 
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,0 
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1 


2 
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1 
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L050N256HISEED 
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,0 
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1 
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1 
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0.001 
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0.001 
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1 


1 
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1 
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0.001 


100 
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1 
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0.15 


0.001 


100 


1 


1 



present in seed-mass BHs (grey curves); Fig. [5] shows the 
redshift zero ttibh — M« and mBH — cr relations. We asso- 
ciate BHs with gravitationally bound objects by identifying 
bound su bstructures in the simulation using the algorithm 
SUBFIND l|Springel et al.ll200ll : iDolag et al.ll2008l ). We note 



that in this plot we show total halo stellar mass as a function 
of BH mass, as opposed to the observations, where only the 
bulge stellar mass is calculated. This means that all curves 
can be shifted slightly to the left. Finally, Fig. |9] shows the 
median specific SFR (SSFR) in bins of stellar mass. In this 



14 C. M. Booth & J. Schaye 





Figure 5. Successively zoomed projections of the gas density in a 3 Mpc/h thick slice from our L100Af512 simulation at redshift zero. 
BHs are represented in this plot by open circles and the area of each circle is proportional to the logarithm of the BH mass. The largest 
circle in the lower-left panel represents a BH of mass 3 X 10^ Mq 



plot, the grey lines represent results from simulations that 
do not include AGN feedback. In figures [9] and |8] the vertical 
lines represent the halo stellar masses at with 50% and 90% 
of haloes contain BHs massive enough to have performed at 
least one heating event. It is immediately clear from Fig. [8] 
that the msH — o" relation is much more robust to changes 
in parameters than the msH — M, relation. 

Each set of simulations is compared to our fiducial sim- 
ulation {L050N256), which uses the model parameters that 
were justified in Sec.|4] To aid comparison between the differ- 
ent simulation sets, the fiducial simulation appears in every 
plot as a solid, black curve. Details of all the simulations 
discussed in this section appear in Table [3] We now discuss 
each simulation set in turn. 



5.2.1 The effect of box size and mass resolution 

We consider first the effect of changing the box size at a con- 
stant resolution by comparing models L100N512, L050N256 
and L025N128. The size of the simulation box has a negligi- 
ble effect on both the star formation rate density (Fig. [6^) 
and, for z < 4, on the global mass of BHs (Fig. [7^). Because 
the properties of individual BHs are set by local physical 
processes, increasing the box size does nothing to the scal- 
ing relations (Fig. [8^) , except for allowing us to probe the 
mass function up to larger halo masses. The same holds true 
for the SSFRs of individual objects (Fig. |9^). 

We now assess the impact of numerical resolution on 
our results. We compare simulations at three different res- 
olutions {L100N512, L100N256, and L WON 128) both with 
and without AGN feedback. Simulation Z/IOOA'^128 has a 
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Figure 6. The dependence of the global star formation history on various sets of parameters, (a): The effect of changing the simu- 
lation volume at a constant resolution; (b):The effect of changing the resolution at a constant box size; ('cj: Simulations with different 
prescriptions for the generation of seed BHs (mhalo.mim "iseed) \ (d)'- Simulations with different values of the heating efficiency {e{);(e): 
The effect of changing the way in which AGN feedback energy is distributed to the surrounding gas particles {AT^i^, rihcat); (f)'- The 
effect of changing the BH accretion model (qq, /3). Curves in grey represent simulations that do not include AGN feedback and the solid 
black curve in each panel represents the fiducial simulation (L050N256). 



dark matter particle mass of 3.54 X 10^° Mq, a factor of 64 
worse than that used in L050N256, our lowest resolution 
production simulation. 

We concentrate first on the star formation history in 
these simulations (Fig. [SJa). At high redshift [z > 4), the 
star formation rate is governed by numerical resolution. At 
low redshift [z < 2) the two highest resolution simulations 
{L100N512 and L100N256) have star formation rates that 
differ by ~ 0.2 dex, indicating that this quantity is not yet 



fully resolved. We see by comparing simulations with and 
without AGN that the factor by which AGN feedback de- 
creases the global SFR is the same in both simulations. How- 
ever, in the lowest resolution simulation {L100N128; blue, 
dashed line), the AGN are largely ineffective at decreasing 
the global SFR. 

We now turn our attention to the BH properties. The 
z < 1 integrated BH density (Fig. [TJj) is virtually indis- 
tinguishable between runs L100N512 and L100N256. This 
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Figure 7. The dependence of the evolution of the the cosmic BH density on the parameters of the AGN model. The panels are the same 
as in Fig. [6l the shaded grey area represents the redshift zero BH density as measured by Shankar et al. (2004). The grey curves show 
the cumulative density in seed BHs. The solid, black curve in each panel represents the fiducial simulation {L050N256). 



occurs because the global black hole density is dominated 
by the massive BHs, which are well resolved in both simula- 
tions. L100N128 underpredicts the redshift zero BH density 
by a factor of two. This is due to two reasons. Firstly, in this 
simulation seed mass BHs are placed only in objects with 
masses larger than 3.54 x 10^^ M©, which means that we ne- 
glect a large population of important black holes. Secondly, 
we see that the first BHs in this simulation start growing 
only at low redshift (2 « 4), and as such, may not have had 
enough time to grow onto the observed BH scaling relations. 

This picture is borne out by the properties of individual 
BHs. Fig. [HJa shows that for a > lOOkm/s the BH proper- 



ties are well converged, whereas at the lower mass end of 
the scaling relations the resolution becomes important. The 
lowest resolution simulation underpredicts the BH scaling 
relations at all masses. 

If we now examine the SSFRs of individual objects 
(Fig. [9]d), looking first at simulations that do not contain 
AGN {L100N512NOAGN and L100N256NOAGN- grey, 
dotted and solid lines respectively) we see that the SSFR 
is almost converged above a stellar mass of lO^^ '^Mo. If we 
now consider the effect of adding AGN feedback to these 
simulations, we see that the stellar mass at which the simu- 
lations with and without AGN diverge from one another has 
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Figure 8. The Magorrian relation and the rviBH — relation for each set of simulations. Line styles and panels correspond to those 
in Fig. |6] The simulated stellar velocity dispersion is calculated as the mean of the three one-dimensional velocity dispersions of the 
stellar particles in each halo. The pale, solid, grey line in each plot represents the redshift zero observations of Haring & Rix (2004) (left 
panel) and Tremaine et al. (2002) (right panel). The solid, black curve in each panel represents the fiducial simulation {L050N256). The 
horizontal, triple-dot-dashed line indicates the mass that a BH requires to begin to heat its surroundings, and the horizontal grey, dashed 
line indicates the black hole seed mass in the fiducial simulation. The two solid vertical lines show, for the fiducial simulation, the stellar 
mass above which 50% and 90% of haloes contain a BH massive enough to have performed at least one heating event. 



not yet converged, indicating that results on the scale of in- 
dividual objects in these simulations may remain affected by 
numerical resolution. We note that the redshift zero global 
SFR is affected by resolution to a lesser degree than the 
SSFR. This is likely to be because the total galaxy stellar 
mass is an integral over all time of the galaxy star forma- 
tion rates, so the strong resolution effects at high redshift 
(Fig. |6)d) persist in the redshift zero population. 



We thus conclude that, because local processes govern 
the size of individual BHs, the simulation box size is unim- 
portant when discussing BH scaling relations. However, the 
limited simulation mass resolution means that the stellar 
masses in our simulated galaxies are not yet fully resolved. 
We also conclude that BH properties, such as the integrated 
cosmic mass density and BH scaling relations are converged 
in all of our production simulations, but decreasing the mass 



18 C. M. Booth & J. Schaye 




Figure 9. The median SSFR as a function of galaxy stellar mass for each set of simulations. Line styles and panels correspond to those in 
Fig. [6] The pale, grey curve in each plot represents the SSFR as measured in the simulation without AGN feedback (L050N256NOAGN). 
The solid black curve in each panel represents the fiducial simulation (L050N256). The upturn in galaxy SSFR in the simulations without 
AGN feedback at M* fa 3 X 10^" Mq is due to SN feedback becoming inefficient. The two vertical lines in each panel indicate, for the 
fiducial simulation, the stellar mass above which 50% and 90% of haloes contain a BH massive enough to have performed at least one 
heating event. 



resolution by a factor of 64 places us in the regime where 
no AGN properties are resolved. Therefore, we can conclude 
that the global results we present in this paper are robust 
to changes in numerical resolution, but that predictions in- 
volving the stellar properties of individual objects should be 
treated with caution. 



5.2.2 The seed model 

Two parameters control the way that seed BHs are gener- 
ated: the BH seed mass {rrisocd) and the halo mass above 
which BH seeds are inserted into haloes (mhaio.min)- For the 
minimum halo mass we chose to use a fixed number of par- 
ticles, mhaio,min ~ lOOmDM , rather than a fixed halo mass 
because this ensures that BHs are always placed into well 
defined haloes, regardless of the numerical resolution. Con- 
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sideration of the AGN feedback process then allowed us to 
place a minimum value on the possible mass of BH seeds 
(10~^mg) such that the average time between heating events 
for a BH accreting at the Eddington rate is shorter than the 
Salpeter time, making it possible for the BH to regulate its 
growth. We also noted that the requirement that the total 
mass in seed mass BHs generated in a simulation is much 
lower than the observed redshift zero cosmic BH density, 
does not allow the seed mass to be much larger than the 
minimum possible value. We thus chose mscod = 10"'^ mg as 
our fiducial value. We now examine how changes in these 
two parameters affect our results. 

By comparing our base model {L050N256) to a simu- 
lation in which BH seeds are only placed into haloes that 
are ten times more massive {L050N256HIHALO; green dot- 
dashed curves) we can say that most of our results are insen- 
sitive to the choice of mhaio.min- The global star formation 
rate density in the simulations (Fig. [6j;) is virtually unaf- 
fected by changing mhaio.min from 10^ ttidm to lO^mDM, 
probably because the SFR of galaxies is only affected by 
AGN in the most massive haloes (M. > lO^Mg; Fig. 
The cosmic BH density is also insensitive to mhaio.min- This 
follows because in the fiducial model, the cumulative seed 
BH mass makes up only ~ 2% of the redshift zero BH den- 
sity, so removing some seeds has a negligible effect on the 
global density. We now turn our attention to the BH scaling 
relations (Fig. |8};) , where we see that the BHs grow on to 
the observed scaling relations at a higher mass, but that for 
M* > 10^^ M© or (T > 100 km/s the results are almost iden- 
tical to the fiducial simulation. These are the BHs that are 
energetically important and we can therefore conclude that 
results derived from our AGN simulations are insensitive to 
uncertainties in the minimum halo mass that contains a BH. 

Large changes in the mass of the seed BHs 
{L050N256HISEED, L050N256LOSEED) can lead to more 
significant changes in the global properties of the simulation. 
If the initial seed mass is lowered by a factor of ten then by 
redshift zero the global BH density is slightly greater than 
in the fiducial case (Fig. [T];). This occurs because initially 
smaller BHs take longer to grow onto the BH scaling re- 
lations both because they need to grow more and because 
they grow more slowly (see Fig. [1} . As a result of this the 
galaxy potential well is deeper by the time that the BH 
begins to heat gas, and so it requires more energy input 
for the BH to stop its own growth. Because it takes longer 
for AGN feedback to become effective, the global SFR is 
higher (Fig. |6j;). The same argument can be made in re- 
verse for high seed masses and explains why increasing the 
BH seed mass slightly decreases the global star formation 
rate (Fig. [6j;). We can draw the same conclusions from ex- 
amining the SSFR of individual objects (Fig.|9};). Changing 
the seed mass has virtually no effect on the SFR of the most 
massive objects, and only affects haloes near the threshold 
mass at which AGN feedback begins to become important. 

The cumulative amount of seed BHs in the simulation 
L050N256HISEED is within a factor of four of the redshift 
zero black hole density. Although this does not violate our 
constraint that the total mass in seed BHs should be less 
than the observed redshift zero value, it is barely satisfied, 
and leaves very little room for AGN in our simulations to 
grow. Additionally a large value for the initial seed mass 
means that we place seed BHs initially above the BH scal- 



ing relations (Fig. [8];) . Nevertheless they still grow onto the 
same scaling relations. 

We note that all our global results are a very weak func- 
tion of the seed mass. For example, changing the seed mass 
by a factor of 100 changes the global SFR by no more than 
a factor 2.5. Other uncertainties, particularly those intro- 
duced by our lack of knowledge about the way in which 
BHs accrete, lead to us conclude that the specifics of the 
seed model are relatively unimportant. 

5.2.3 The feedback efficiency 

The efficiency with which a BH's radiated energy is coupled 
to the ambient gas, ef, is treated as a free parameter in our 
simulations. We show in this section that this parameter 
controls the normalization of the total mass in BHs an d of 
the BH scaling relations (see also lDi Matteo et al.ll200"5l) . In 
our simulations ef was tuned to match the z = ttibh— »Tihaio 
relation and the redshift zero cosmic black hole density. A 
value of Ef = 0.15 was found to work well. 

Changing tf from its fiducial value has no discernible 
effect on the global star formation rate (Fig. [6ji), or on 
the SSFRs of individual objects (Fig. |9}1), but the total 
redshift zero BH density is inversely proportional to the 
feedback efficiency (Fig. [7}1). Most strikingly, simulations 
L050N256VHIEPS, L050N256HIEPS, L050N256LOEPS, 
and L050N256VLOEPS have values of et that differ by fac- 
tors of 4 , 2, 1/2, and 1/4 from the fiducial run, and the ratio 
of their z = Q BH densities to that of the fiducial simulation 
are 0.24, 0.53, 1.93 and 4.03 respectively. The fact that the 
final mass of BHs is directly proportional to 1/ef implies that 
in our models any given BH grows until it has injected a spe- 
cific amount of energy per unit halo mass, at which point it is 
able to reduce its local density a nd to effectively s e lf-reg ulate 
its growth. As demonstrated in lBooth fc Schayg (120091 ) this 
remarkable r esult ag r ees q ualit atively with t he id eas pre- 
sented in e.g. ' Fabianl (|l999i ') and lSilk fc Ree s' (1998), where 
BHs grow until they can expel gas from the galaxy, at which 
point they enter a quiescent phase. 

The redshift zero global BH density (Fig. [7}i) and the 
normalization of the BH scaling relations (Fig. |8ji) are both 
directly proportional to 1/ef, and so we can use these obser- 
vations to constrain the value of ef that allows our model to 
be consistent with redshift zero observations. We employ a 
value of ef = 0.15 in the present study, but note that we will 
show elsewhere that the parameters of the AGN model must 
be tuned in conjunction with the SN feedback prescription. 

5.2.4 The heating mechanism 

The resolution of cosmological simulations is too poor to re- 
solve the scales on which AGN inject energy into the ISM 
and as such we are forced to make some assumptions. It is 
therefore important to verify that the results obtained from 
our simulations do not depend strongly on the implementa- 
tion of the energy coupling mechanism. 

Two parameters control how energy is distributed from 
the AGN to its surroundings: the minimum temperature in- 
crease and the number of neighbouring SPH particles to 
heat (ATmin = 10* K and ihcat = 1)- The fiducial value 
of the minimum temperature increase was constrained by 
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two competing effects. Firstly, if ATmin is too low, the cool- 
ing time of heated gas will be so short that the gas will be 
able to radiate away its energy before having a dynamical 
effect. On the other hand, if ATmin is set too high, then the 
amount of energy needed to perform a heating event will 
be so high that the numerically determined time between 
heating events exceeds the Salpeter time-scale, making it 
impossible for the BH to regulate its growth. We found that 
ATmin = 10* K provides a good compromise between these 
two effects. Similarly, to minimize the numerically controlled 
duty cycle we set riheat = 1 in our fiducial models. Here we 
assess the impact of distributing feedback energy in different 
ways. 

If we change either the temperature to which gas is 
heated from 10* K to 10^ K {L050N256T7) or the number of 
black hole neighbours affected by each feedback event from 
1 to 10 {L050N256HINHEAT), we see small changes in the 
cosmic star formation history (Fig. |6^) and in the cosmic 
BH density (Fig. [7^). We note, however, that the results are 
relatively insensitive to these parameters, and that changing 
the value of either ATmin or riheat by an order of magnitude 
affects both BH properties and SFRs by only ~ 0.3 dex. We 
thus conclude that as long as the feedback efficiency (ef) is 
calibrated such that the redshift zero black hole relations 
are satisfied, and that the minimum heating temperature is 
sufficient for feedback to be effective, our other results are 
robust to the precise way in which this energy is coupled to 
the ISM. 

5.2.5 The accretion model 

Because numerical simulations cannot resolve the properties 
of the multi-phase ISM, we can justify the use of high values 
for the multiplicative factor, a, in the Bondi-Hoyle accretion 
rate (see Eq. [TJ, at least for high density gas (Sec.|4]). This 
provides a motivation for the constant-a accretion model - 
similar to those previously published - where BHs accrete 
with a very high efficiency. In this paper we have introduced 
a new class of models that treat accretion of low-density 
gas differently. By noting that the ISM will not contain a 
cold (T <C 10* K) phase at low pressures and that we are 
also able to resolve Bondi-Hoyle accretion onto BHs with 
masses greater than the resolution limit of our simulation if 
we resolve the Jeans mass, we argue that a must asymptote 
to unity in low density environments. We then parametrize 
our ignorance about the state of the multiphase ISM, and the 
precise mechanism by which AGN accrete by introducing a 
parameter, /3, that describes a power-law dependence of the 
accretion rate on the local gas density (Eq.|4]). 

The growth of a black hole depends upon the accre- 
tion model used. In the constant-a model (Bondi-Hoyle with 
ao = 100) the growth is Eddington limited unless the gas 
in the immediate surroundings of the BH has a density that 
is much lower than typical of the ISM, e.g. as a result of 
feedback from the BH. In the constant-/? models, however, 
much larger densities are required for the accretion rate to 
become Eddington limited. Efficient growth is only possible 
if the BH's local density is enhanced by dynamical processes. 
BHs can only decrease their accretion rate (and hence reg- 
ulate their own growth) when densities are low enough that 
accretion rates are no longer Eddington limited. The den- 
sity at which BH accretion rates become Eddington limited 



depends on the accretion model (Fig. [2]). In the constant-a 
model with ao = 100 BHs self-regulate at densities below the 
density at which a multi-phase medium is expected to form 
(Fig. ^ , this has a large effect on the physical properties of 
the galaxy and because the simulations resolve Bondi-Hoyle 
accretion in this regime it invalidates the assumptions used 
to justify the use of ao = 100 in the first place. 

We now consider how changes in the accretion model 
affect the properties of the simulated galaxy population. In- 
creasing P from 2 to 4 (simulation L050N256B4) depresses 
the global star formation rate somewhat (Fig. [Bf) because 
the j3 parameter controls the gas density at which the accre- 
tion rate for a given BH mass becomes Eddington limited 
(Fig. [2}. Hence, for a given value of /? there is a critical value 
of the local gas density above which a BH can 'switch on' 
and begin to grow rapidly. For larger values of l3 this occurs 
at a lower density, and hence in smaller haloes. This mani- 
fests itself in the BH scaling relations by changing the mini- 
mum galaxy mass at which the BHs grow onto the observed 
scaling relations. A lower value of /3 therefore increases the 
global star formation rate by allowing BHs to grow only in 
more massive haloes. Increasing (3 above 4 does not have a 
large effect on any of our results, as for any large value of 
P, the BH accretion model behaves in such a way as to be 
Eddington limited in star- forming gas (Fig. O, and the dif- 
ference in behaviour between a /3 = 4 model and a /3 = oo 
model is very small. Any physical process that is strongly 
dependent on the local density is affected by numerical res- 
olution, because higher resolution simulations allow the for- 
mation of higher density regions. For this reason we caution 
that the stellar mass at which BHs grow onto the observed 
relation is affected by numerical resolution, and so the value 
of l3 may need to be tuned to different values for simula- 
tions with mass resolutions significantly different to those 
presented here. 

A comparison of the effect of different (3 models on the 
SSFR of haloes (Fig. [Sf) supports this picture. The stellar 
mass above which the behaviour of each AGN model di- 
verges from the behaviour of the simulation without AGN 
depends upon the value of (3. 

In the constant-Q model with ao — 100 
{L050N256A100B0) BHs grow in an Eddington lim- 
ited manner from their birth, and as such suppress star 
formation in all haloes that contain a BH. This is visible in 
Fig. [of, in which the constant-a simulation deviates from 
the simulation without AGN in haloes with a low stellar 
mass. The constant-a simulations efficiently suppress star 
formation in every halo that contains a BH and as such the 
integrated SFR is more than an order of magnitude lower 
than in the other simulations. The same effect is present 
in the large simulation volume, but is less pronounced 
because seed mass BHs are placed only into larger haloes, 
where both accretion models are capable of suppressing 
SF. This result implies that in order for a simulation that 
employs a constant-a accretion model to reproduce the 
observed break in galax y properties at log(M,) ^ 10.5 (e.g. 
iKauffmann et al]|2003l ). the resolution must be tuned such 
that BH seeds are placed into haloes of the correct siz^f]. 



^ Because the accretion rate also depends on the assumed effec- 
tive EOS, this may not be true for all constant-Q models used in 
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Considering now reducing the value of ao in the 
constant-a models from 100 to 1 {L050N256A100B0 com- 
pared to L050N256B0; a constant-/? model with /3 = is 
equivalent to a constant-a model with ao = 1), equivalent 
to removing the numerical 'fudge factor' from the BH ac- 
cretion rates, we see that BHs are unable to grow in all but 
the very most massive haloes (Fig. |8f), and so the halo SS- 
FRs are virtually the same as for the simulation without 
AGN feedback in all haloes up to a stellar mass of 10^^ Mq 
(Fig. 13). This, in turn, means that the global SFR density 
in the ao = 1 simulation is very similar to that in the simu- 
lation run without AGN (Fig. [3). The global BH density at 
redshift zero is actually higher than that in the fiducial sim- 
ulation (Fig. (Tf) ■ This occurs because the BHs that do grow 
in the ao = 1 simulation grow very late, and so find them- 
selves in a deeper potential well. In order to self-regulate 
they then need to output more energy and hence grow even 
more than in the other simulations. The global BH density 
is therefore by itself not a good indicator of how efficiently 
AGN are able to suppress star formation. 

We now consider increasing ao by an order of magnitude 
relative to the usual value {L050N256A1000B0 compared 
to L050N256A100B0). We already know that the density 
below which BHs can self regulate depends strongly on the 
value of ao (Fig. [2J, and it follows that the value of ao 
will have a strong effect on the physical conditions in the 
centres of haloes. In the case of ao = 1000 the global BH 
density is similar to the ao — 100 case (Fig. [7f), whereas 
the global SFR and galaxy SSFRs are much lower than in 
any other simulation (Fig. [Sf and Fig. [9f), and the mBH — 
Tistoiiar relation is shifted significantly to the left as the BHs 
effectively shut down star formation in all haloes. However, 
the mBH — o" relation is not as strongly affected by the very 
effective AGN feedback. This is likely due to the fact that 
the stellar velocity dispersion tracks the depth of the galaxy 
potential well, which contains a large contribution from dark 
matter, whereas the galaxy stellar mass is sensitive to the 
details of the feedback processes operating inside the galaxy. 
We see that changing the value of the free parameter ao can 
have profound effects on the simulated galaxy population, 
even if the evolution in the global mass density of BHs is 
barely affected by the parameter change. 

We conclude this section by noting that predictions 
from AGN feedback models of this type are sensitive to 
the accretion model that is used. More generally, it is 
not clear that the Bondi-Hoyle rate is the correct accre- 
tion rate to use in the case o f a BH accreting from a 
hot halo (|Krumholz et al.| [2006;) . We find that a different 
parametrization of the accretion rate leads to profound dif- 
ferences in the star formation histories and speed of BH 
growth and therefore caution the reader that the AGN ac- 
cretion mechanism represents a significant source of uncer- 
tainty in all our results. 



the literature. In the models of 305, which employ an EOS that is 
initially stifTer than our 7^3 = 4/3, the BH growth time has a lo- 
cal minimum at njj = which suppresses AGN growth relative 
to our constant-a model (Springel, private communication). 



6 DISCUSSION AND CONCLUSIONS 

We have presented and tested a method to incorporate 
SMBHs into cosmological, smoothed particle hydrodynam- 
ics simulations. The method, which is a substantially mod- 
ified version of the one introduced by SOS, self-consistently 
describes the mass growth of BHs and feedback from AGN. 
Here we consider growth through mergers with other BHs 
as well as through accretion of gas. The AGN feedback in 
our model is thermal and local to the BH. 

Although we also use the SPH code GADGET III, our 
code differs from that of S05 in many ways, including the 
use of different models for star formation, feedback from SN, 
radiative cooling and stellar evolution. Particularly relevant 
for AGN feedback is the fact that, contrary to SOS, we do not 
make use of a subgrid model to describe the different phases 
of the ISM. Following SOS, we make use of subgrid BHs to 
allow BH masses to be small compared with the masses of 
the particles containing the BHs. We note, however, that 
while this approach allows one to significantly extend the 
range of BH masses in the simulation, the accretion radius 
will only be resolved if the BH mass exceeds the local Jeans 
mass. Unfortunately, this is generally not the case if the BH 
mass is small compared with the particle mass. 

The AGN model is specified by seven main parameters 
(Table [l|. Two parameters describe the BH seed generation 
mechanism: the BH seed mass, mgeod, and the halo mass 
into which seed BHs are placed, mhaio.min. Two parameters 
describe the amount of energy that is coupled back to the 
ISM per unit accreted mass: tr, the radiative efficiency of 
a BH accretion disk, is the fraction of the rest mass energy 
accreted by the BH that is radiated by the AGN and tf 
is the fraction of the radiated energy that is coupled ther- 
mally to the ISM. For a given BH accretion rate, the rate at 
which energy is injected into the ISM depends only on the 
product Ertf. However, we do need to specify the radiative 
efficiency since the mass growth of the BH is proportional 
to (1 — Er). A further two parameters control the numerical 
implementation of the injection of energy into the ISM by 
AGN: the number of neighbouring gas particles heated by 
each BH heating event, nheat, and the minimum amount by 
which their temperature is increased, ATmin. Because we let 
BHs store feedback energy until they have saved enough to 
heat rihoat particles by ATmin degrees Kelvin, these last two 
parameters together determine the AGN duty cycle for a 
fixed accretion rate. Finally, we require one additional pa- 
rameter that controls how BHs accrete gas, and we describe 
two different models for this process. 

The gas accretion rate is assumed to scale as the Bondi- 
Hoyle accretion rate, evaluated at the location of the BH and 
on the scales resolved by the simulation. We do not, however, 
allow the accretion rate to exceed the Eddington rate. The 
Bondi-Hoyle accretion rate predicted by the simulation will 
underestimate the true rate if the density is underestimated 
or if the temperature is overestimated. A lack of numerical 
resolution may result in an underestimate of the gas den- 
sity, which motivated SOS to multiply the Bondi-Hoyle rate 
predicted by the simulation by a constant factor ao = 100. 
We call models of this type constant-a models. We noted 
that cosmological simulations lack not only the resolution, 
but also the physics to model the cold ISM phase. For ex- 
ample, our simulations use a poly tropic effective equation of 
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state for densities at which the gas is expected to be multi- 
phase. Hence, they will miss the cold component for which 
the Bondi-Hoyle accretion rate would be highest. This will 
lead us to strongly underestimate the Bondi-Hoyle rate in 
high density gas. 

Although cosmological simulations cannot yet resolve 
cold, interstellar gas, many do resolve the Jeans scales at 
densities low enough for the ambient ultraviolet radiation 
to suppress cooling below 10'* K. Specifically, any simula- 
tion that resolves the Jeans scales in the gas surrounding a 
BH particle, will also resolve the Bondi-Hoyle radius if the 
BH mass exceeds the local Jeans mass. Hence, at sufficiently 
low densities the Bondi-Hoyle accretion rate is modeled cor- 
rectly and multiplying it by ao ~ 100 would result in a large 
overestimate of the accretion rate. 

We therefore introduced a second class of models in 
which the Bondi-Hoyle accretion rate is multiplied by a fac- 
tor (nn/nn)'' for densities nn > n^, where njj is the density 
above which the gas is expected to be multiphase (we take 
nn =0.1 cm~^). We refer to this class of models as constant- 
/? models. Note that both constant-a and constant-/? models 
use a single free parameter. Because we have changed the 
density-dependence of the accretion rate, we cannot claim 
to be simulating Bondi-Hoyle accretion, even though the 
changes are motivated by the Bondi-Hoyle formula and even 
though we do retain the Bondi-Hoyle scaling with the BH 
mass. We argued, however, that this is also true for constant- 
a models because the use of values ao 3> 1 implies that the 
densities and temperatures predicted by the simulations are 
greatly in error. 

The parameters ao and /3, used in the constant-o? and 
constant-/? models respectively, control the ambient gas den- 
sity at which the BH accretion rate becomes Eddington lim- 
ited (Fig. [2}. Because the maximum densities sampled by 
the simulation increase with halo mass (both because more 
massive haloes are resolved with more particles and because 
the central pressure increases with the depth of the poten- 
tial), these parameters effectively set the halo mass above 
which BHs can grow efficiently. We set P — 2, which we find 
results in efficient BH growth in haloes with stellar masses 
> 10*°'® Mq in these simulations. Using a constant-Q pre- 
scription with ao = 100 implies that, in the absence of AGN 
feedback, the accretion rate is essentially always Eddington 
limited. Because the accretion is efficient even at relatively 
low gas densities, AGN feedback is in that case important in 
all haloes that exceed mhaio.min. For this class of models the 
halo mass above which AGN can suppress star formation 
is thus in effect a free parameter (mhaio,min). For BHs with 
masses greater than 10^ Mq self-regulation can only occur 
at densities orders of magnitude below the star formation 
threshold (Fig. In this regime we resolve Bondi-Hoyle 
accretion, invalidating the assumption used to justify large 
values of a in the first place. Constant-a models therefore 
underestimate the gas density required for self-regulation 
and will thus overestimate the suppression of the star for- 
mation rate. 

Having chosen a prescription for gas accretion, we then 
derive values for the other model parameters. Because each 
BH stores its feedback energy until it suffices to heat ^heat 
neighbours by ATmin, we are faced with a numerically de- 
termined duty cycle (for a given accretion rate) . In order for 
the BH to be able to regulate its growth, we require the time 



between heating events, Atheat oc riheat ATmin, to be as small 
as possible and certainly smaller than the Salpeter time if 
the accretion is Eddington limited. We use ATmin ~ 10* K, 
which ensures the temperature of the heated gas is high 
enough that the injected thermal energy is not just radi- 
ated away, and riheat = 1, which minimizes Atheat. Because 
Afheat decreases with the ratio of the mass of the BH to 
that of the heated gas particle, the requirement that Athcat 
is smaller than the Salpeter time for Eddington-limited ac- 
cretion implies that the (subgrid) BH mass must exceed 0.1 
per cent of the gas particle mass (Fig. [3]). Hence, we set 
rrisced ~ 10~^mg. We set rrihaio.min = lOOmoM in order 
to ensure that seed BHs are placed only into well defined 
haloes. These parameter values will obviously need to be in- 
creased if they result in seed and/or minimum halo masses 
that are lower than expected physically, as may be the case 
for simulations that use a much higher resolution than is 
used here. 

We assume the standard value tr = 0.1 for the radiative 
efficiency and tune et to match the redshift zero msH — JWhaio 
relation and cosmic BH density. A value of ef = 0.15 was 
found to provide a good match to the observations. 

Having specified the AGN model, we then analysed the 
results from a large suite of cosmological simulations chosen 
to investigate the sensitivity of the predictions on the model 
parameters. For this purpose we compared the predictions 
for the cosmic SF history, the cosmic BH density, the red- 
shift zero BH scaling relations (both the Mbh — M* and the 
Mbh — o" relations), and the redshift zero galaxy specific star 
formation rates (SSFRs). 

We demonstrated that the fiducial model provides good 
agreement with both the observed mass density in BHs 
(Fig. [7]) and the BH scaling relations (Fig. |8} , and that the 
inclusion of AGN feedback in the simulations effectively sup- 
presses star formation in galaxies with stellar masses greater 
than > 10^°'^ Modot (Fig. (9)1. We will discuss the comparison 
between the simulated global SFR density and observations 
elsewhere, but for now we note that the SN feedback pa- 
rameters of our models were tuned such that a simulation 
without AGN feedback has a peak SFR density that is in 
good agreement with observations. As such, adding an extra 
source of feedback energy inevitably results in an underesti- 
mate of the SFR density. In order to achieve a good match 
with observations the properties of the SN model must be 
tuned in conjunction with those of the AGN model. 

However, the focus of our study was not to match ob- 
servations, but to explore the dependence of the results on 
the parameters of the model. Our main conclusions are sum- 
marised below: 

• Regardless of whether BH seeds are initially placed 
above or below the BH scaling relations, they grow onto 
the same relations. 

• Because the global BH density is dominated by massive 
BHs and because AGN feedback is only important in high- 
mass haloes, uncertainties in the seed model employed do 
not lead to significant changes to the global properties of our 
simulations. Changing the initial seed mass by two orders of 
magnitude changes the global star formation rate by only 
a factor of 2.5. The assumed seed generation model can, 
however, affect galaxy properties at around the galaxy mass 
where AGN first begin to become energetically important. 
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At higher masses galaxy properties are largely insensitive to 
the initial distribution of BH seeds. 

• A s discussed more comprehensively in lBooth fc Schav3 

the normalization of both the global BH density and 
the BH scaling relations is nearly exactly inversely propor- 
tional to the AGN feedback efficiency, ertf- Most strikingly, 
changing the efficiency by a factor 16 does not give rise 
to any discernible changes in the global SF history. These 
results imply that the total amount of thermal energy in- 
jected by AGN is conserved when the feedback efficiency is 
changed. These results can be explained if BHs grow until 
they have generated enough energy to regulate the accretion 
rate and if the required amount of energy depends on the 
depth of the potential well on scales that are larger than the 
radius on which the BH dominates. 

• Changing the way in which the thermal energy from 
the AGN is distributed in the gas surrounding the BH has 
little effect on our results, as long as the gas is heated to 
a temperature that is high enough for its cooling time to 
become long. Hence, thermal feedback can be efficient in 
cosmological simulations that do not resolve the multiphase 
ISM. 

• Cosmological simulations currently underestimate the 
Bondi-Hoyle accretion rate in dense gas because they lack 
both the resolution and the physics to model the dense, cold 
phase of the ISM. It is therefore necessary to increase the 
predicted accretion rate by a fudge factor, either by explic- 
itly multiplying the accretion rate by a numerical correc- 
tion factor or by employing a subgrid model for the unre- 
solved gas physics to artificially boost accretion rates in star- 
forming gas. Using a multiplicative factor that asymptotes 
to unity in the regime where the simulation is able to model 
the relevant physics (our constant-/? model) gives different 
results from using a high factor throughout (the constant-a 
model) as has been used in most previous work. In general, 
the density above which BHs are able to accrete efficiently 
depends upon the accretion model used (Fig. Because 
higher mass haloes contain higher density gas, the accretion 
model determines the halo mass above which AGN feedback 
becomes effective. Until the simulations are able to resolve 
Bondi-Hoyle accretion in a multiphase ISM, the predictions 
of the models will remain subject to significant uncertainty. 

• The TTiBH — o" relation is more robust than the ttibh — 
relation to changes in the model parameters (Fig. [8]), and so 
provides a test of the numerical model that is less affected 
by uncertainties in numerical parameters than the Magor- 
rian relation. This is likely due to the fact that the stellar 
velocity dispersion tracks the depth of the galaxy potential 
well, which contains a large contribution from dark matter, 
whereas the galaxy stellar mass is sensitive to the details of 
the feedback processes operating inside of the galaxy. 



In summary, we have presented and tested a new model 
for the self-consistent growth of BHs and feedback from 
AGN in cosmological simulations. In a future work we will 
discuss the interaction between AGN feedback and other 
physical processes, and show that the results obtained from 
an AGN model depend also on other processes such as SN 
feedback and stellar mass loss. 
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